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1.0 INTRODUCTION 

This document reports on the development of two new shell finite elements 
for application to large deflection problems. The elements in question are 
doubly curved and of triangular and quadrilateral planform. They are re- 
stricted to small strains of elastic materials, and can accommodate large 
rotations . 

A principal goal of this work has been to obtain accuracy for strongly non- 
linear problems commensurate with that which is obtainable in linear finite 
element work. This particular goal presents a special challenge, because 
the presence of strong nonlinearities in plate and shell finite elements 
of conventional formulation cases excessively stiff behavior. 

The elements described in this document make use of a new displacement func- 
tion approach specifically designed for strongly nonlinear problems. In 
order to accommodate this approach, two further new developments were nec- 
essary, relative to strain displacement equations and nonlinear solution pro- 
cedures. Thus, the developments reported herein represent an attempt at a 
considerable extension of the current state of the art for nonlinear two 
dimensional elements. 

The two doubly curved shell elements developed, a triangle and a quadri- 
lateral, are based on the relatively simple linear elements of References 2 
and 3. The displacement function development for nonlinear applications is 
based on the beam element formulations of Reference 1. The procedure is 
called herein the HMN method, in recognition of its first use by the authors 
of Reference 1, The strain-displacement equations developed are of the shallow 
shell type specially tailored to the HMN procedure. Additional terms have 
been included in these equations in an attempt to avoid the large errors 
characteristic of shallow shell elements in certain types of problems (e.g., 
Ref. 7). An incremental nonlinear solution procedure specifically adopted 
to the element formulation was developed. This solution procedure is of 
combined Incremental and total Lagrangian type, and uses a new updating 
scheme. A computer program was written to evaluate the developed formulations. 
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This program can accommodate small element groups in arbitrary arrangements. 
Two simple problems have been successfully solved, and are reported in Sec- 
tion 10. The results indicate that this new type of element has definite 
promise and should be a fruitful area for further research. 



2.0 GEOMETRY AND COORDINATE SYSTEMS 


There are three coordinate systems that are required for the two new elements; the 
local coordinate system, the global coordinate system, and the solution coordinate 
system. The elemental matrices are derived in the local coordinate system and are 
eventually merged and solved in the solution coordinate system. With slight varia- 
tions, these matrices can be discussed for both the triangular and quadrilateral elements. 
The geometry of the elements is input in the global coordinate system. The triangular 
element is assumed to be a constant curvature element, whereas the initial geometry 
of the quadrilateral element has the shape functions of the assumed displacements. 

Base planes for the local coordinate systems are established as shown in Figure 1 . 

The triangle has the x axis parallel to the line foining nodes 1 and 2 with positive 

sense from 1 to 2. The y axis is perpendicular to the x axis and lies in the plane 

of nodes i, 2 and 3 with positive sense in the direction of node 3. The z axis is 
determined by the right hand rule. 

For the quadrilateral, a mean plane is established through the midpoints of the straight 
lines joining the corner nodes. The x axis is parallel to the line bisecting sides 
3,5 and 1,7 with positive sense from side 1,7 to side 3,5. The y axis is perpendicu- 
lar to the x axis, lies in the mean plane, with positive sense from side 1,3 to side 

7,5. The z axis is determined by the right hand rule. 

With these definitions, the transformation matrices can be written from the local to 
the global coordinate system. These transformations are presented in the Appendix, 
Section 1 2.1 . . The transformations that are required, however, are from the local to 
the solution coordinate system. The solution coordinate system is selected to minimize 
potential numerical problems by solving the equations in a coordinate system in which 
the z axis is nearly normal to the shell midsurface. This is accomplished by averaging 
the normal unit vectors of the local coordinate systems of the elements adjoining each 
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node. Figure 2. The x and y axes are selected to lie in planes parallel to 
the XZ and YZ planes, respectively. With these assumptions, the transformation 
matrix from the local to the solution coordinate system can be written 

O*) 

where the elements ofj^A] are given in the Appendix, Section 1 2 . 1 . 


[a] M 
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Figure 2: SOLUTION COORDINATE SYSTEM (SHOWN FOR 
TRIANGULAR ELEMENTS ONLY) 
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3.0 SHELL STRAIN - DISPLACEMENT RELATIONSHIPS 


Shell strain displacement equations take a wide variety of forms, depending on 
their intended use, coordinate reference systems employed, and specific terms dropped 
in accordance with thinness or shallowness assumptions. Experience has shown that 
it is permissible to take considerable liberties in the formulation of these equations, 
in order to achieve simplicity, so long as the omitted terms are truly negligible for 
the problem types under consideration. It is important in this regard to perform the 
mathematical derivation carefully, and to make approximations in the final results, 
where the full implications of the approximations are visible. This is the approach 
taken here. For the present work a formulation applicable to stepwise nonlinear 
analysis of a shell element is required. This application suggests several specific goals 
for the strain formulation: (1) To handle large rotation effects, the transverse shears 
should contribute to in -surface, as well as norma I -to -surface equilibrium; (2) To avoid 
undue complexity in the strain equations due to accumulating shell midsurface deforma- 
tions and associated nonlinearities, element local coordinate systems and incremental 
displacements should be cartesian; (3) To facilitate implementation of an HMN type 
constraint procedure, each increment of displacement should involve small rotations, 
and updating of local cartesian coordinates should be carried out. These goals can be 
achieved by deriving exact strain displacement equations of a shallow shell element in 
cartesian coordinates, and dropping appropriate terms in the final result. The approach 
taken is a combination of the incremental and total Lagrangian approaches since in- 
cremental Lagrangian strain increments are used, and in addition total Lagrangian strain 
is used in equilibrium checks during the nonlinear solution procedure. 

Kirchhoff Shell Theory 

Consider the shell element shown in Figure 3 . The element is referred to two 
coordinate systems. The x-y-z system is cartesian, and the x-y plane is called the 
local baseplane of the element. The element in its initial unstrained state is 
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displaced an amount CJ° from the x-y plane. The projection of the x-y system 

fK f\ 

onto the midsurface of the shell element in its initial position defines the X ,^j 
system, according to the equations; 

x°= X 

y; 3 . C3-i) 


r\ 

where the superscripts denote the initial state of the element. 2 is perpendicular 

A ^ A 

to X , (j . It is noted that X and ^ are slightly distorted scales of true length. 

a 

Points of the shell element located a distance 2 above the midsurface of the shell 
have the x-y-z locations 


0 n 
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A 

0/ 

X 21 x 

- 2 


" X 

-z a 

a ^ 

A 


A 


0 

~ 2 

A 

* y 

- zco 


Z = 6J' i-Z 


( 3 - 1 ) 


After a deflection u, v, w, the coordinates of points on the shell are given by 

0 • 

X = X+a-2 6J 

Lj -- </ +-U-5 oj' C 3 ’" 3 ) 

The directions of u, v, w are, respectively, the positive x, y, z axes. Thus, 
cartesian displacements rather than the customary ones parallel and normal to the shell 
surface are employed. 

Before deformation, the base vectors of the curvilinear x, y, z system are 


with 
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_ — , etc. 

5? " 0? 


( 3 - 4 -) 




(3-s; 
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I0P1?0 




where i, j, k are the unit cartesian base vectors, 
and, using the above notational shorthand. 


J-Sv"), (cj'V 

3; '[(S.cj"'') , (I -£ cj°"), (cj°‘)J O 6 ) 

3 i - ( 1 )] 


The quantities set off by commas are, respectively, the x, y, and z components of 
the vectors. After deformation the base vectors are 


°\ C = 


~Si + [(u'-Zoi"), (y ’-iej’'); & 'J - g'j + 

3S ~ ?J +L ( “ ) i (*' -£&"), ts'J ~ *■ ^ 

9a = 9s (Wj, ifo)7 = 9s "T s 


M 


The metric tensor of the X ,Z coordinate system is defined by the scalar products 
of the base vectors. 


$5 = $ ; * 3° 

(s-s) 

3'J ~ 3; ' Sj ~ (?i + %;)' (jj ♦ A J 

J 

(?-3) 
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Here, the subscripts denote X , ij , or Z . It can be shown that for small deforma- 
tions the strains which occur in passing from the undeformed to the deformed state 
are described by (Ref. 4 ) 


Thus, 


. a 




which becomes 


£ ,j = * [% • 7 


Performing the indicated operations yields the results for the strain tensor for a 
Kirchhoff shell theory 


, y X 

a + a 


a') & + (yO e ^ (cj ') 2 J ^ 

- Z [cj- + u ' (<xr°"+ oj') + v (oj- 0,/ i~ <*J ') 

-if 7 


Cqq ~ + * ^-T°' Cjjr’ 

* J 

-z[ur"-r v f (cj°"*-cj“) + u (cj° '‘+<J m/ )J 

-z*[ 7 • 
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2 cfj j = 2 <5 
x y 9* 


V +6t -f U U * V y + CJ CJ +- \CJ- Cif + GJ CJ j 

-7 ’ ( * %. o// ) -t. ' /, L. s r°*f ^ 

— J / a / r / *f* / i \ /. T /. T / 

£ u uy ***- <-«- v My / *-*- j w-' / 

+-V‘(ej'' , -f'CJ c ' / ) + v'( 

-r[ 7 


These equations are to be used incrementally, and the cartesian x, y, z system 
for each increment will be obtained by updating immediately prior to the increment*. 
Under these conditions, the displacements u and v will be very small compared 
to w, and the squares of their derivatives can be "neglected. The terms which are 
neglectable are underlined. However, in the bending terms, the terms involving 
products of u and v derivatives with w and w° derivatives will be retained, 
as these terms provide x- and y-direction forces due to transverse shears, the absence 
of which in shallow shell theory can cause serious errors. With the approximations 
indicated, the incremental equations are 



s A CL (cf J + AJ ) A 

- z[acj'+ Aa (<J°"gj") + U (a Cj") + &Y ((J + & 

+ v'(a(J 0/ )J 




This updating is not a change in element shape; details are discussed later 
in Section 8. 
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A V (uj° \ cJ ') A CJ* " 

H [a cj“ + & V'( CJ°"+ca") j ' v'Acj" +■ Au'((J° m/ + kj ') 
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(3-18) 


In these equations, w denotes the displacement normal to the xy plane, accum- 
ulated prior to the incremental displacement under consideration. Note that 
while the direct strains are actually physical strains, the shear strain 
tensor values are half of the physical strain. 

Consider Figure 4. Here a relatively large displacement of a beam element is 
shown. The element is approximately inextensional, so that accompanying the 
large w are necessarily large u displacements. It can be shown easily that 
3u/3x is also large. If u and w are referred to the updated base plane which 
corresponds to the end of the displacement and denoted, respectively by u and 
w, then u becomes small and 3u/3x becomes negligible. This reasoning justifies 
the use of the above strain equations, with the underlined terms neglected, 
for the calculation of total strains corresponding to the end of an increment 
and the associated updated base plane. This is important to avoid cumulative 
updating of strains, which can cause relatively large errors. 
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Shell With Transverse Shear Deformation 

Figure 5 illustrates for an x-z section view the deformations considered in this 
section. Again, u and w are cartesian. The rotation about the y axis, 
controls the x i -direcrion displacements of material points above and below the shell 
midsurface, according to the right hand rule convention. The coordinates of points 
on the shell are given Ijy 2 y , ay° , since no ^ or 0 y need be associ- 
ated with the initial shape, w , After deformation, the material point coordinates are 


The base vectors are 


X =- X + U + Z 


9 = 


UJ + % 


9? ” 9x ), "2^ ) , ' ] 

9y"99 + l (u Q)i 

h +[ 0 y 


(s-\d) 


(s-'Zo) 


The strains are 


u'+i. f (u'f * (vT + (uj f] * cj'cj"' 

-i [-a, * u (u“-A,') * V '(»“'*&) 

' -it ] 

- v‘ + 2. J (ll ) T (v’) Z ^ ( U) ) J ^ OJ cj ° 
-z[ei+tr(*j°--e;) * v(cj°'-+e;)] 

. .. J 




(j-'Z'l) 
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« - C ^ - 2 [a + v '1 +■ z[ Uu + v\y + oj oj'J ^ a [cj cj + cj <uj j 

*< 4 yx L J t- 


_ i~ +&* * u '( ) t u - / (^°’ e tj) 

+-v'(cj°' f +e;)^v'(tj^9;)] 

2 Z £ . 1 

- i i r (" uj0 + Q + ^ - ^)" • ' • * 

^ *- I; 9 y)' h %. ('CJ°-@*) " J 


£ gz = 2 r + 



(V*4) 


As before, the underlined terms and the higher order terms in Z are neglected. 

In this non-Kirchhoffian case there exists a transverse 
shear strain whose contribution to the energy of deformation will have to be handled 
carefully (Ref. 5 ), because of its tendency to cause excessive stiffness. 

The incremental equations are, dropping the noted terms 


A6gg 


= Au'+ (cj°'+ cj') A 6J z[-AB s + Au' (v°"- s *) * Av '< 



- u'A@y ■* v'a @*J 



= Au(0J°"-&i) + 4V‘ 

(ur + S, ) 


-a'A0, + V'A<%] 
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~ a. A O x 


As discussed above for the Kirchoffian case, total Lagrangian strains can be computed 
from the above , provided that updating of the element baseplane coordinate system, 
with the corresponding transformation of the total displacement components, is carried out. 

Shallowness and Baseplane Derivatives 

/\ A ^ 

Recalling the definition of the coordinate system X , Q , 2 , it is seen that the 

r\ r\ r\ ^ 

mapping of the x-y system into the X - system has caused X and X, and also 

/\ 

and Lj , to have identical numerical values. Thus, all differentiations with respect 
to X and in the strain equations have identically the same values if they are re- 
placed by, respectively, x and y differentiations. This is required because element 
displacement functions will be expressed as functions of x and y. 

In addition, for the strain tensor values given to be physical strains, it is necessary 
that the metric tensor components and g^ be essentially unity. That is, it 

is necessary for the lengths in the shell surface to be approximately unity; 


-/ , i +* (cj**y 


^3 -3l) 
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Thus, in identifying physical values with the strain tensor components, it is 
necessary that 

(ar°T> 1 

For the work reported herein, this assumption is made based on the fact that iach 

•T# 

element, referred to its local baseplane, is very shallow, both before and after 
deformation. Thus, all strain tensor values given are physical values. Note, 
however, that the physical shear strains are twice the tensor values, . 


£ 
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4.0 VIRTUAL WORK FORMULATION 


The formulation of the elemental stiffness matrix which includes the HMN degrees 
of freedom and the degrees of freedom to be reduced by minimum energy considera- 
tions, is virtually identical for the triangular and quadrilateral shell elements. The 
general formulation, proceeding from the virtual work expression, is presented in this 
section. The details of individual matrices required by the particular element are 
presented in Sections 6,0 and 7.0. 

The virtual work is given by 

where the strain is expressed in terms of the displacement gradients as 

✓ 

or 

fe; = (RO;‘ +-R 1 ;j k G K ) ST 9j 


( 4 - 1 ) 

( 4 -' 2 -) 
( 4 - 3 ) 


For the triangle, the displacements can be written directly in terms of the generalized 
degrees of freedom. The quadrilateral element, however, requires an additional 
Jacobian type transformation from the displacement gradients in the parent element to the 
displacement gradients in the local coordinate system. In general this can be written 
as 


Bj -GO 




GOj„ - GOj„ 


for the triangle 
for the quadrilateral 


(4-4) 


Thus, the strains can be written 

- §,v 


(4-5) 
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B ir = (ffO;j *Rl;' k 6> k )G0j, 

(4-c) 

- (ro,; * R h k doi^, * ) 

The stresses are 

or ^a’*c>o u Ae j , 

£ ' 

where O' is the initial stress and £)<9,/ is the matrix of material constants. 

Substituting th e above in the virtual work equation yields 

) S s : dH ‘ l 4 ~ Q ) 


The incremental load deflection equation is required 


AP;'- K0 ;jX 

KO;; -j' [ 4 ; D0 si +ao m X-%„A Po k , 5 14-10) 

where higher order products of the incremental displacements have been ignored. 

This gross stiffness matrix includes the deformational degrees of freedom, , plus excess de- 
grees of freedom associated with the minimum energy procedure, , and degrees of 
freedom associated with the HMN constraints, oQ, Changing from indicia! notation 
to matrix notation, the generalized degrees of freedom can be written in partitioned 
form as 

(4*11) 

The triangular element is derived with several options as to the number of degrees 
of freedom to be constrained by HMN constraints. Consequently, the generalized 
displacements cx^ assume several forms to account for this variability. This is 
discussed in greater detail in Section 6,2. 
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The stiffness equation then has the form 


’ ll<o] Uf] 


C4-12.) 


In general, the HMN constraint can be written as 

, r-v. -i T /i C* *>. 

The HMN constrained stiffness matrix is then 

m-i?] r Mtc] 


(4-13.) 


L4 - 1 +) 


Finally, the minimum energy degrees of freedom , are reduced out by standard 

procedures to yield the stiffness matrix 

This matrix relates the generalized degrees of freedom to the associated generalized 
forces. For the quadrilateral the generalized degrees of freedom a can be written 
directly in terms of the nodal degrees of freedom, ^ . The triangle, however, Is 
derived in terms of the deformational degrees of freedom and these degrees of freedom 
must be supplemented by the rigid body degrees of freedom 

? * [T*] ff] ^-iO 

where m -Ex 3 for the quadrilateral. 

Applying this transformation yields the elemental stiffness matrix 

[K ] -[T*l r [K*Uri ivn) 

relating elemental nodal forces to elemental nodal degrees of freedom. 
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5.0 HMN Methodology 


5.1 Background 

The incorporation of nonlinear terms in the strain displacement equdtions can give rise 
to problems of excessive stiffness of finite element models for some large deflection 
problems. This difficulty was discussed first by Mallett(Ref 6 ), and later by 

Haftka, Mallett, and Nachbar (Ref 1 ). These authors found that poor problem 
solution accuracy was obtained for large deflection and post buckling problems of beams, 
and formulated a procedure for eliminating the problem. They called the resulting 
element a "stability element." 


The formulation of Haftka, Mallett and Nachbar (HMN method ) can be simply 
explained as follows. The axial strain for small rotations and strains for a beam element 
is given by the formula. 


c , i / 'I 2, 

c* = 57 *T(sr) l s-V) 

Customarily, u is linear in x, and Uf is cubic. The functional form of £ K is therefore: 
from u, constant; from UT , through the fourth degree In x. This situation holds also for 
the linearized incremental case. 




<^4u 

Jx 


Jur c)A ur 

dx 3* 


(s-0 


through the product 


c)oJ 

Jx 


\nA6j) 
1\~7TJ * 


The result of this situation is that, if the £ x strains due to r— are retained in their 


then 


d* 


1 II VJ I t 

entirety , /superimposed on the constant strain due to ~ are in general, first, second, 

third, and fourth degree-in-x terms due to The strain energy, since it is a 

Jx 
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quadratic form, is necessarily increased by these additional terms, even if their average 
should be zero. This fact, of course, clearly indicates that the element, and associated 
finite element models, are too stiff. 


The formulation given by Haftka, Mallett, and Nachbar avoids this difficulty by intro- 
ducing supplemental u displacements which include functions of the second through the 
fifth degree in x. Through the term , these functions create strains of the first 

through fourth degrees, which are sufficient to compensate the element for the higher 
degree \-j^ J type terms. The authors impose the axial equilibrium equation t 


JM* 

5k 




(*- 3 ) 


to solve for the amplitudes of the supplemental a functions. This use of the equilibrium 
equation is equivalent to minimizing the total potential energy of the element, i.e., 
to an elemental level reduction of the supplemental u freedoms. 

It is noted that in the case of the beam element no inter-element incompatibilities can 
result from the above equilibrium, or minimum energy, conditions. This results from 
the fact that the supplemental functions cause no nodal displacement, i.e., the 
supplemental functions involve only internal displacements. In attempting to apply the 
HMN procedure to plate or shelf elements, this latter convenience does not occur. For 
these cases, the high degree nonlinear strains due to the OT derivatives occur not only 
within the element, but also on its edges. Thus the supplemental U and V' functions 
must be nonzero both inside the element and on its edges, and the matter of inter- 
element compatibility must be considered. In addition, the high degree nonlinear 
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strains which are to be "controlled" by the supplemental functions are now functions of 
two spatial coordinates instead of only one, with considerable attendant difficulty. 

For these reasons, the details of the application of the HMN procedure to plate and 
shell elements are not straightforward. The overall procedure is reasonably clearcut, 
however, and is briefly explained below. 


5.2 Selection of HMN Functions 


We began by illustrating the selection of HMN functions by a simple example, and 
then proceed to a more general discussion. In cartesian coordinates, the nonlinear 
OJ derivative terms are fir) ' (jfi ' and (&r)( sir) ' contributing, respectively, 
to , 'Ey , and • Similar terms apply to the incremental strains, with the 

linearized term being, for example, These types of terms create 

high degree polynomial strains illustrated by the following several equations, which 
are representative for a finite element bending formulation 

CJ= (i-y)(}-x z ') 

(jr- * 4 -) 






(tt)** ( i -3.J +<]')(**) — * £. 


x nonlinear 


In this example, the nonlinear is of the second degree in tj and X 
that the basic element a. displacement function has the form 

UL- ( oc-t- ) (<L + dt ) 


Presuming 


(*■-*•) 
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The basic strain form is 


* b(c.*ty) 


Ls-t) 


Thus the only strain which can be provided by the basic • function is the one 
linear in (j above, and the supplemental u. function, here called a , must provide 
the terms listed below if it is to cancel out completely the higher degree nonlinear 
strains. 

~ 3 $23 

U. ~X , J X 


(f-7) 


c 

The u functions are of course one degree higher in X than those of 

^ 3 

These u. functions are seen to contain the lowest degree in X member, X . 

This is a consequence of the example chosen for illustration, in which the oJ function 

also 

is not complete. In an actual element derivation, 6/ would/bontain also linear in X 

terms, and ’ ’> these in combination with the second degree term would result 

^ a 

strain terms containing X and hence the necessity for ti terms like X . 

This simple example suggests a procedure for choosing the polynominal forms of the 

|»"V» 

supplemental u and v functions. The steps involved are: 

1. Identify polynomial forms of basic element strain £< , , 

2. Identify polynomial forms of nonlinear strains due to uf , A OT 

3. Construct supplemental U, , v such that the combinations CL+ u j v + V 
contain strains of polynomial degrees through the highest nonlinear 
forms created by UT , A ur . 
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4. Examine the constructed functions to assure that they form a complete 
set with no intermediate omitted terms, and augment the function set 
if necessary to assure this. 

This last step poses particular difficulties, which are discussed briefly below. 




The procedure above can be followed rather easily for any type of element under 
development. In general, however, it must be followed by an additional task which is 
highly element dependent, and here lies the essential difficulty in HMN function 
development. The difficulty lies in the fact that u and V functions must be constructed 
which, in combination with u. and v , form a complete set; and which also do not con- 
tribute to the displacements, such as nodal freedoms, which characterize the basic 
element U. and V functions. In general this means that U and v must consist partly 
of purely infernal (pillow) functions, plus other forms which accomplish the needed 
deformations of the sides of the element, without causing any changes in the values of 

, r%t 

the basic element freedoms. The pillow functions are needed to make the total U-+-LC 

/Nf A/ /V* 

and Vfy sets complete. The values of U and V along the sides are also subject to 
the completeness requirement, as well as the requirement to be independent by U and 
V . Thus the w and V forms, which as developed in the procedure given are simple 
polynomial forms, must be manipulated info element dependent forms suitable for the 
particular element under development. This is best illustrated by two examples, given 
here for functions of a single variable, but capable of generalization to the case of 
two spatial variables. Figure 6 shows a case in which the basic U. function is 
cubic and a must contain terms through the fifth degree, due to a cubic form of „ 

This ?s a case in which the element freedoms are nodal displacements and first 
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derivatives. The task here is to construct fourth and fifth degree functions which do 
not affect values of the displacements and their first derivatives at the ends x = + 1 . 

They have been constructed simply by inspection, and are shown in the figure. Note 
that all polynominal degrees from the zeroth through the fifth are represented by the 

rv 

combined u and U, , in order to obtain a complete set of functions. 

Figure 7 is a case in which the element basic freedoms make use of a midside node. 
The basic a and CJ functions are quadratic, and a cubic u. is required. For elements 
belonging to a family with nodes located along the sides, for example the isoparametric 
family, it is most convenient to construct the higher degree supplemental functions by 
conceptually adding more nodes along the sides. Thus, in Figure 7 , the basic 

element has a single midside node and quadratic displacement functions. The addition 
of the supplemental is to be accomplished by adding more side nodes. A cubic U. 
is needed, and hence two nodes in the interior of the side would be sufficient. This 
is not a desirable arrangement, however, because the center node should be retained. 

The addition of two nodes proves workable. The resulting set of five nodes on a side, 
including the corner nodes, is capable of representing a fourth degree polynominal, 
whereas only a cubic is required. This is not a real difficulty, however, and the 
indicated cubic shape for u. provides the necessary function. This shape has the 
necessary zero displacement of the end and center nodes. Thus, in using additional 
side nodes as an aid to construct the supplementary functions, it is not necessary to 
form the highest degree function of which the nodal arrangement is capable. 


Suppose as a further example that a quartic were needed. In this case the two quartic 
diagrams shown in the figure would be used. These two shapes together contain a cubic 
and a symmetrical quartic, as is easily verified by forming the quantities 

- *(x~4:Vx-£) ^ 

The use of these functions therefore makes unnecessary the explicit use of the cubic form. 
This example shows the advantage of the procedure of constructing HMN functions by 
adding midside nodes consistent with the nodal arrangements of an isoparametric family. 
This procedure has readily provided cubic and quartic forms, while simultaneously 
omitting all terms of the basic displacement functions, simply by assigning displacements 
to all nodes other than those which participate in the basic function shapes. 

The matter of the required pillow functions is usually resolvable based on consideration 
of what pillow functions would be required for a higher order finite element whose basic 

rv /v* 

displacement functions include polynominal forms of degree comparable to cl and V . 
Thus for the second case above, pillow functions corresponding to a higher order (5 nodes 
per side) isoparametric element would be used. For the first case above, if for example 
a triangular element based on area coordinc+es is being developed, the known pillow 
functions of higher degree area coordinate forms can be used. 

The above discussion has illustrated by example two methods for constructing HMN 
supplementary functions. One method is applicable to cases in which element freedoms 
include displacements and displacement derivatives at corner nodes, and the other two 
cases in which displacements of both corner and intra-side nodes are used as basic 



freedoms. In each case, the procedure begins by constructing functions of the needed 
polynomial degree along the sides of the element, while causing no displacements of 
the basic element type, and achieves completeness of the LM-U and V V seteby adding 
pillow functions as needed. This procedure meets the conditions which, at the present 
stage of development, appear necessary in HMN function methodology. 
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5.3 HMN Constraints 


The added displacement freedoms introduced as HMN functions must be constrained or 
solved for as a part of the finite element formulation or solution process. There are 
several alternatives here, and generally the approach to this problem is not straight- 
forward. Some of the relevant details are discussed in this section, and further detail 
is given in the later sections dealing with the two shell elements developed in the subject 
research. 

The use of the HMN functions is a recognition of the need for highly competent membrane 
displacement functions in nonlinear bending problems. These functions can be thought 
of as simply added displacement freedoms, and subjected to inter-element compatibility 
requirements and solved for in the global solution process. In this case the HMN con- 
straint methodology would not exist as such. However, this approach has the drawback 
of increasing problem size; it is desirable to keep to a minimum the number of global 
freedoms, in addition, the philosophy underlying the HMN functions — that these 
functions are to remove complex forms and keep the strains in the element simple — 
suggests that the HMN functions be dealt with on the elemental level and that the basic 
element membrane functions are sufficient in the global solution process. 


31 


REPRODUCIBILITY OF THE 
ORIGINAL PAGE IS POOL 


Elemental level HMN constraints can be of two types; minimum total potential energy 
(equilibrium); and restrictions on the functional forms of the membrane strains. The two 
approaches can be used together in a single element derivation, operating on different 
HMN freedoms. In both cases, consideration must be given to the implied inter-element 
displacement incompatibilities which are implied by the constraint procedure. The 
specific strain constraint- is considered first. 


Suppose that the membrane strains, say <5^ , due to, respectively, the basic element Li 
functions, the HMN U functions, and the bending oJ functions are given by 
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In these forms, x and y are spatial variables and the double subscripted U. / U- , and 
CJ are coefficients which depend on the element freedoms of the type U., u , and 
Oj , respectively. The upper limits of the sums have been set by the choices of basic 
element functions and the HMN functions. Note that for the first summation, the 
limits M.| and N^ are small numbers tike 1 or 2, while the values and Ng are 
larger, like 4 or 5, and that these upper limit values are the same for the second and 
third summations. The latter condition it will be recalled was set by the choice of the 
HMN functions. 


One approach to the HMN constraint procedure is to simply equate coefficients 
and C4*„ for those strain contributions which ore of higher degree than M^ and N^ , 
i .e. , 


^ m n ~ ^yn/\ \ W ^ 

(5 '-lo) 

U* y y )/) _ b-fyn y\ j K7 ^ 


In the first case above, while M r*| may be any value; and in the second case the 

similar condition on m and n occurs. This approach generates a large number of 
equations and generally involves duplicate constraints, i.e., singularities of the con- 
straint equations. It has been tried and found unworkable in the cases attempted. The 
procedure amounts to a complete elimination of high degree strains everywhere in the 
element (above M^ and N^ levels). It appears that such an approach would be success- 
ful, and it is not clear why it failed to work out. It is recommended that this approach 
be given consideration in possible future work on HMN elements, since alternative approaches 
provide only a partial elimination of the high degree strains, and thus would appear to be 
inferior. 
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Another, simpler approach il to eliminate the high degree strains along specific lines in 
the elements. For example, the high degree yC forms can be eliminated along a set of 
lines e * c ‘ This takes the form 


t.(CL, )x%E = 03 

n-o n-o ^ 

n*o n-o ^ 


(5*1.1) 


This results in fewer and more easily formulated constraint equations and is a much simpler 
procedure overall than the one first described. If enough lines, i.e,, values y^, y^ t 
etc. , are used, it should give an excellent control of the high degree strains. In order 
to relate the two procedures described above, it is noted that if the strain due to U. and 
CJ , say , is known to vary as a quadratic in tj , then the constraint should be 

% 

duce singularities in the constraint equations. 

This example involves control of the % functional forms of the x-direction strain <kx; 
the lines in the element are in the same direction as the strain itself. Of course, 

r .... . . f\. 

additionally could also be controlled in the behavior tj along specific lines 
X~ X', j Xe. 3 etc. But if the number of iji lines chosen is adequate to represent the 
known functional behavior, as illustrated above for a quadratic form, then this two- 
direction control of would likely (or certainly) produce singularities in the constraint 
equations . 


, and ij should be used, 
j values necessarily pro- 


applied along three tj = constant lines, i.e"., <j - <■ 
Fewer U> values fail to give complete control, and 
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For the control of the shear strain, e , a similar procedure can be used, 
and in this case the two-directional control has been used without intro- 
ducing singularities. Thus, for the polynomial forms y 11 are controlled 

along x=x^,X2 . •• , and the forms x n are controlled along the lines y=y^,y2 
This is done for the values of n which exceed those provided by simply the 
basic element membrane displacement functions. 


Control of follows the analogous pattern to that of e^. 

The above described procedures refer to rectangular elements, for which the 

directions of the strains and the directions of the element sides coincide. 

For triangular elements a procedure has been used in which the HMN functions 

control strains along the sides of the element and are again handled as 

described above. Hence, if s denotes the coordinate along a side and n the 

coordinate normal to a side, the strain e and £ can be controlled along 

s sn 

all three sides. Isoparametric quadrilateral elements also deal with strain 
directions which differ from the directions of the element sides. In this 
case a feasible procedure appears to be to impose the HMN constraints on the 
"parent" element, and carry the derived constraint equations on for use in 
the subsequent distorted elements. The above completes the discussion of 
HMN function constraints based on explicit consideration of the forms of 
the strain functions. The procedures discussed are illustrative in nature 
rather than steps of an established methodology. Inventiveness has been 
found to be required for two dimensional HMN element development. 

The imposing of HMN function constraints by means of minimum energy 
condition is a formal methodology which can be applied without particular 
difficulty. It is briefly discussed below. This constraint amounts to a 

set of equilibrium 
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conditions, and is therefore the type of approach originally adopted for beams by the 
authors of inferences 1 and 6 . However, for two dimensional elements it 

raises questions of inter-element displacement compatibility. These are discussed in 
Section 5.4. 


The minimum energy control of the high degree strains amounts to an elemental level 

/v _ 

reduction of the U. and v freedoms. Thus, if the complete element stiffness matrix is 
partitioned in the form 


k = 


fe* J k 
h, ! k 


(S-VL) 


u J \ 
I J 


where the lower rows pertain to the U and V freedoms, the constrained stiffness 
matrix refering to the basic elements u , V , and OJ functions above has the familiar 
form 

Du - /?;: v] (s-13) 

This approach tends to produce an overly flexible element due to inter-element displace- 
ment incompatibilities. 

5.4 Inter-Element Displacement Compatibility 

HMN function- constraintseither by explicit strain function control or by the minimum 
energy conditions may create inter-element displacement incompatibilities. This 
difficulty can be largely avoided in the former procedure by a careful selection of the 
specific strain controls imposed, at the expense of relaxing some of the desired controls 
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on the high degree strdin functions. It cannot be so directly avoided in the minimum 
energy approach. This section discusses the inter-eiement compatibility matter and 
provides recommendations for preferred strain control constraints. 


Consider Figure 8 . Line a-b is the side of an element, and coordinate axes x and y are 

perpendicular and parallel to a-b.- The basic element functions are denoted by u, v, and 
W, and have amplitudes defined by the nodal freedoms. Consider the determination of 
the supplementary function V , which involves displacements parallel to the side a-b. 

The strain is completely controlled by the displacements V , V and the slope pp 
The value of V determined by HMN constraints imposed on is completely controlled 

(aJ~ 

by pp , since the HMN constraints only apply to those high degree polynomial j 

forms which are unaffected by the lower order basic function \/ . Hence the v function 

along a-b will be identical for the elements on the left and on the right of a-b, because 

-jp is necessarily the same for these two elements. Still further information concerning 

V is available, however. Since V is determined not only along a-b, but also throughout 

the adjoining elements, for example by constraints imposed on along lines 3 a r W 

o>i7 

in the figure, it is possible to draw conclusions concerning . Not only does the 


correspondence 


C5--14-) 


hold, but, provided V has been given x-direction variation comparable to that of p7 
also the first term Taylor expansion holds, and 


, ’■[^(sr)] 'X (5--15) 
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Thus, H can be concluded that 


/dv\ 

^ / 

f <)(*! ) 

( ** Lt, 

m 

, I 


.-U 


O- 


in which the order of differentiation has been exchanged. Since is 

continuous across a-b, ^ is also continuous, it can be concluded that -j~ 




<5v/ 


as well as V is continuous across the line a-b. !t is noted that can be discon- 
tinuous across a-b; this is required to allow equilibrium for the case of a y-direcfion 

dv 

line load along a-b. The continuity of across a-b will in general be approximate, 
as it depends on the competence of V as a function of x in each element and also on 
the exactness with which continuity of ^ is enforced. Nevertheless, clearly 


tends toward continuity between adjoining elements, across x - constant lines. 


Consider next HMN constraints on the shear strain • The relevant high degree forms 

contain contributions from ) an< ^ ^ rc?m an< ^ * Continuity between 


th 


. dcj 


e elements on the left and on the right of line a-b occurs exactly fot -j — and approxi- 
mately or exactly for , depending on the particular finite element formulation 

employed. Hence — + is essentially continuous across a-b. It has been seen 
^ ck 

that -Tfir is at least approximately continuous, and hence also is continuous 

across a-b. Since — is continuous across a-b, and since U vanishes at the nodes a and 
b, U itself is necessarily continuous across a-b, or at !ea e * approximately so* 
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It has been seen that HMN constrai ntson and C iL ^ along x = constant lines pro- 
vides exact or approximate compatibility of Lk and V between elements across these 
same lines as inter-element boundaries. It remains to consider the use of an HMN con- 
straint on ;„ x along an x = constant line, for example line a-b in Figure 8 . This 

strain constraint relates and SSL . Since the latter is continuous, or cpproxi- 

mately so, across a-b, this constraint makes approximately continuous across 

a-b. This particular result is not really essential, inter-element continuity of LA and 
V already having been obtained, and hence the use of such an HMN constraint, i.e., 
on an extensional strain whose fiber direction is normal to the line along which the con- 
straint is imposed, is not recommended. The likelihood of dependent constraint equations 
together with the lack of basic need for this constraint are the basis for this 
recommendation . 


One further matter concerning explicit strain constraints and inter-element compatibility 
needs discussion. Consider again the constraint on along line a-b in Figure 8 

- / 'l f c ^’ r ) 

For the higher degree polynomial: forms in y, the strain is defined by ( / 


and -—+■ / an d the discussion given earlier applies. Howevenfor the hiphest 

a<j ©u " 

degree x basic function in e and the lowest degree x HMN function in e , 

xy xy 


due, respectively, to the terms 




and 


d x 


, there is a special 


problem. For this particular poiynominal degree only / an HMN constraint condition 
would involve the correspondence 


c)u. <5^ 
*r — t — 


d Uj\f aar'j 


dX 
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in which both basic functions and HMN functions are involved. The term sliL is not 

dy 

in general continuous across line a-b in the figure, and hence for this particular 
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polynomial term neither would — be continuous. The major difficulty in 
application of this constraint is that it is inconsistent with the constraint based on <Sy 
and involving in a form which is a function of x. As discussed earlier, the 

constraint tends to match up the x-directional functional forms of and F 

while imposing equality of their y-direction polynomial coefficients. The lowest 
degree HMN constraint terms oppose this x-functional form match-up, and 

thus can be expected to cause difficulties in the HMN constraint equation set. It 
is recommended that the lowest degree HMN constraint be omitted. No inter- 

element displacement compatibilities will be caused by this omission. 


Finally, the matter of constraint of the HMN functions by minimum energy conditions 
requires discussion. In this procedure, the element is effectively made as flexible as 
possible by setting the w* and v freedoms to values which minimize the strain energy. 
It is virtually guaranteed that the result of this constraint applied to HMN edge 
functions will be inter-element displacement discontinuities. The possibility is strong 
in this case that an excessively flexible element will result. Therefore, minimum 
energy constraints on HMN functions are only recommended for the internal (pillow 
function) freedoms of the HMN set. 

The above rationale for the inter-element continuity of V due to an C Sn constraint 
is dependent on the continuity of ■— . If all the HMN constraints are as des- 
cribed, with, in particular, conditions imposed on lines parallel to the 
side , this continuity should be obtained. However, it may occur that the internal, 
pillow-type HMN functions are controlled not by explicit strain constraints but by 

* Here v is normal to the side, u is parallel to the side, consistent with 
definition used in the element derivation. 
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minimum energy conditions. In this case, if these pillow functions are capable of 

du 

non -zero on the element edges, then the reasoning leading to continuity 

/V/ _ 

of v based on an constraint becomes less certain. The difficulty is that 

constraints are operating on ^ which do not necessarily produce inter- 

element continuity. It would appear that even in this case, if explicit con- 
straints are imposed along lines parallel to the sides, then continuity of and 

hence of v , through an £g n constraint, would be obtained. For the triangular 
element described in Section 6, however, this is not the case. Here, constraints are only 
used along the sides themselves, ^ anc | minimum energy conditions control 
the pillow functions. In this case it would have to be determined by numerical 
experimentation whether the combination of £ s and C SYl constraints, with minimum 
energy internal function constraints, yields satisfactory results. In the initial formu- 
lation of the triangle done in this contract, not only 6 S and £ Sy) but also 8 n 
explicit strain constraints were imposed along the sides. It has been verified that this 
was unsatisfactory. Inter-element incompatibilities of v occurred, and the resulting 
element was too flexible. Consequently, V was omitted along the sides of the element. 

5.5 Total HMN Functions and Constraints 

The foregoing has tacitly dealt with application of the HMN procedures to incremental 
finite element calculations. A subsequent section of this report will deal with iterative 
corrections of incremental solutions in order to avoid accumulating errors due to step- 
wise linearization. It is noted here that the HMN function magnitudes are also subject 
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to accumulating error, and that they also should be improved through iterative 
corrections. All the HMN methodology given in the earlier sections applies equally 
to total or incremental strains, and thus is applicable to iterative calculations. 

5.6 Stiffness Matrix Development Sequence 

HMN constraints can be applied at the outset of the elemental stiffness matrix calculation, 
or they can be applied after the elemental matrix is calculated. In the first case, the 
stiffness matrix never involves, explicitly, the HMN freedoms, as they have been made 
dependent on the basic element freedoms before the matrix generation. In the second 
case, the elemental matrix is initially derived in terms of all the (A , v , UX , U. and 
V freedoms, as though they were independent, and the constraint is accomplished later 
in a standard stiffness matrix transformation. 

Those HMN freedoms which are to be eliminated by minimum energy constraints are always 
handled in a way similar to the second procedure. The elemental matrix is developed 
retaining these freedoms, and their elimination is then accomplished by matrix partitioning 
and a standard elemental level reduction. The second type of procedure has been used in 
the vork of this contract for all HMN function determination. 

The HMN equations and their use to constrain the stiffness matrix are outlined 
briefly below. In the discussion, the unknowns {Act } are the coefficients of 

ci 

the HMN functions used to obtain explicit control of the strains along the 

sides of the element. The unknowns {Aot, } are the coefficients of HMN functions 

b 

which are entirely interior to the element, i.e., which vanish on the sides of 
the element. Since the derivatives of the functions corresponding to the 
{Aa^} do not necessarily vanish on the sides of the element, the specific HMN 
equations involve T Act^ } as well as {AcO. The {Aq*} and {q*} contribute the 
strains due to the bending rotations. 
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In general, the HMN equations can be put in the form 

[«]£*.]-- [ho] J + [Hl(f )] {& % '/ 

where [Hi is a square matrix. Consequently, 


[A «3 = [nr , [no]$^]>inr[nu'}‘-)) Utf 

Expand these matrices by adding columns of zeros such that 





( 4*0 


(9-18) 


Cs-\$) 


(s-^o) 



(s-'ll) 


Then expand with the identity matrix to obtain the constraint equation 



(5 -a'2.) 


which is used to obtain the HMN constrained stiffness matrix. 


[h } -[cVUllc] 
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6.0 TRIANGULAR SHELL ELEMENT 


The goal in developing the triangular element was to obtain a relatively simple 
nonlinear shell element incorporating an HMN capability. The HMN methodology 
requires that relatively high degree membrane displacement functions be used, 
particularly if the basic element normal displacement functions contain high degree 
forms. For this reason, a relatively simple basic element was desired, and the 
BC1Z element was decided upon (Ref 2 ). Simplifications related to a shallow 
shell theory but not subject to all of the shallow shell approximations were desired. 

This was accomplished through the use of a particularly advantageous set of strain 
displacement equations, described earlier in Section 3. 

The HMN methodology was incomplete at the start of this work, having been previously 
formulated only for beam elements. For this reason, a number of options were re- 
tained in the element formulation, relative to HMN functions and constraints. Cal- 
culations have been made to evaluate several of these options. The basic approach 
to the HMN procedure is outlined in Section 5. 


6.1 


ASSUMED DISPLACEMENT SHAPES 


The transverse displacement is assumed to be that of the BClZ element. Reference 
2 . The displacement is composed of two parts, a rigid body translation, w, 

and a deformational displacement w*. 
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The functions F\ have the characteristic that they are zero at each of the nodes 
and the derivatives with respect to x and y are zero at each of the nodes 
except at a node where or Bujx; is nonzero. 


Along each side of the element w is a complete cubic, taking the form of 
the conventional beam bending displacement shape expressed in terms of area 
coordinates. The interior cubic forms are an incomplete cubic. This results 
from a constraint imposed in Reference 2 to assure that constant curvature 
states are available. 
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The membrane type displacements u and v have the same form as the transverse 

displacement w with some additions. The basic constant strain state u, v is 

added to a deformational state u*, v*. In addition, there are interior fifth order 

pillow functions u , v and fourth and fifth order edge functions u , v where 
r m m 

m denotes an edge. Referring to Figure 9 , the displacements are written 
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1<C ej L<, < <J 

^ j i are arbitrary constants 
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The HMN functions are required to control fourth degree strains, and hence 
need to be of the fifth degree in the area coordinates. Along the sides, the 
needed forms are 


r 


and 

(X: o-o 


where ry 1 , y\ take the values 1, 2, 3 and . These functions have zero 

values and zero derivatives at the nodes, and thus do not interfere with the basic 
element freedoms. Each function also is non-zero on only a single side. For 
example, ( "3 , 3 e f~ and (3,d.) . (3 - 3 t ) are non-zero on only the side opposite 
the number 3 node. 


The interior functions required to complete the displacement set through the fifth 
degree are 



(C~T) 


and 

in which and the indices take the values 1,2, 3. These functions and 

their first derivatives are zero at all the nodes. 


The interior functions describe membrane displacements in the element local coordin- 
ate directions. However, the edge functions are used to define displacements referred 
to coordinate axes which are local to the sides themselves. Thus, S and n denote, 
respectively, local side coordinates along and normal to a side, and Up and 
denote the HMN edge functions parallel and normal to the f> th side. The side 
functions and ( 3L 3*)' * (3^-3,,) are used to create four functions; two each 

for each of and Vp , where the subscript of denotes the side opposite 


node IP , and -joivnfn. These edge HMN functions have non-zero values along 
a single side and decrease to zero at the other sides. Each side in turn has 
such functions, aligned parallel and normal to the side. Clearly, no generality 
is lost in choosing the side-oriented directions; they are a matter of convenience in 
implementation of the explicit HMN constraints. 

There are a total of 24 HMN functions, of which 12 are interior, and 12 are side 
functions. Note all of these functions are always used in calculations, however. 

The computer program has various options on HMN function retention and constraints. 
In particular, the side functions V are in some cases omitted. 
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6.2 ELEMENT DERIVATION 


The basic triangular element is derived with several options for handling the excess 
generalized degrees of freedom. Both minimum energy constraints and HMN constraints 
are applied to the excess degrees of freedom to reduce the gross 51 by 51 stiffness 
matrix to the final 27 by 27 elemental stiffness matrix. The program options designate 
which, if any, of the excess degrees of freedom will be HMN constrained, with the 
remaining degrees of freedom reduced through minimum energy. There are presently 
.four, such options available. This will be clarified as the derivation proceeds. 

The virtual work expressions presented ( in Section 4.0 are applicable to both the tri- 
angular and quadrilateral shell elements. The stiffness matrix was obtained in a 
general form in Equation^ which is repeated here for clarity. 




+- 


«... 



(C-9) 


Particular application of this equation to the triangular element displacement functions, 
Section 6.1, and the details of the derivation are contained in this section. 


From Equation (4-4) f the GO matrix has the form 


* [go] [y 7 

L L U x. ; U tj t ) V'y , •> ^ 

Lq - J s LS*;o<*:^J 


S 


*i 




Jr 


v- 


^ osXj ! 


T ar l t 


9 ^i V, O v<jl CJ X Q ( 

^ u > x z. 


U 3 e ux 3 

t = L< J ; Or* L < : or- L^C, l ^2 z °\ 3 J 

L J ~L j , or L J o r L NULL J 


(6- to) 


11 .) 


5Z 


(G- 1 0 


L , J = L ^ a.*- q 6 4 i 
L £ J ~ L cl lt d ,2. o/ 2 , o/g^ d /3 dgj a., A, Q.g a 3 b 3 J 
L °( 3 J = i C-h c 2l C /2 C 22 C /3 A 23 J 


The form of the GO matrix is given in the Appendix, Section 12.2. Similarly, 
the AO and A1 matrices for the triangle can be directly written from Equation^-?.) as 

AE,- {_ftO;j+fll ; . k 0 K ) A 0; Cs-is) 

where 

L AE J "" L A E X x AE yy A £*y A K^ x A Kyy A K jiy J (j^_ 4) 

and L<9J was defined in Equation Explicit forms for AO and A1 are given 

in the Appendix, Section 12.3.. The curvatures (K kkJ Kyy , Kt.y ) are related to the 
bending strains by = Z e.~hc~. Also, physical strains are used instead of tensor strains. 


Finally, the stress -strain relation is given in Equation(4-7)as 

CT;* 00; A A Ej. 


C.G-1S) 


where 

lcr\-iN x A/ y N f y A). M, M, y J 
L &°1 = LN x ' A/ y * At; J 


(G-1g) 


and the matrix of material constants [DO] is given in the Appendix, Section 12,4, 
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With these matrices, the gross stiffness matrix of Equation(c-9)ccin be generated. 

This matrix is of 51 by 51 order involving the deformational degrees of freedom [fj 
and the generalized degrees of freedom [cxCj, \ ck^J. At this point the concept of 
HMN type constraints is introduced. The HMN constraints are applied to the 
degrees of freedom and result in equations of the form 

M = [to] [ff ] * U&j] ^T) 


where 


j I I ® oJ%S- j Q U.T<JZ I 1 Q <J- r j3 


(6-1 8) 


As indicated in Section 5.6, these equations can be put in the form 



-fc] 


'4 ST 

A c< 


and the HMN constrained stiffness matrix is then Equation (fe -7jo) 

m--[?] T MLS] 


(6-13) 


U-'Zo) 


Four options are available for the constraint procedure as indicated in Equations (6-ll) 
and (C-'ii), All of the and [_ degrees of freedom are eliminated either 

through HMN constraints or minimum energy reduction. 


Option 0: 


L0C0.J = Lnull J 

I X 0 

LcS-kJ - L*,: °(rj 

I A 2 V 


No HMN constraints 




Option 1 


L J ~ /. J 

l'°<\ j - 

/ A it 


Only normal strain constraints enforced 


(6 -'2/L) 


Option 2: 


L°< CL J 

/ x 16 

L 

1 x C 


L 0(2 i °^3 J Normal, tangent, and shear strain 
1 , constraints enforced. 

— l_ oC , J 
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Option 3: (j o( J J Only normal strain constraints enforced 

= 1«,| O ; cl, fa, h z a, i>J (G '^ 

|X6 


l<xj 

I x 18 


IX b 


°-Z ^2 C*-’ 
/ x 6 


For Option 3, the six generalized degrees of freedom corresponding to the normal 

e. 

edge displacement, V , are set equal to zero, i.e„, d^, d^/ ^22' ^31 ' 

^32 are zero * The remaining set of displacement functions £L , u T ) u. <2 J v 

•*- X 

V ,v when used in conjunction with the € s HMN constraint yields a fully 
compatible element. The gross stiffness matrix is of order 45 by 45 when the null 
rows and columns have been deleted. 
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Options 0,1 and 3 are subsets of Option 2 ; consequently the most general 
constraints will be developed and the subsets will be obvious. 


The constraint equations are obtained by the procedure developed in Section 5.0. 
As applied to the triangle, the incremental membrane strains along the three edges 
can be written 




Sir? Jpxn-cosoU A I) +[37^*0 


BAoJ- 

In™ 


c Jur 


JAlT 1 


(6-2 S) 


ASL, - 


zn 


1 _ o>AU(ri) 

as 

a cj-° d A or a or 3 Aar 

+- JJTM) -r J^mr 

• 



(«- 


v*» a A u __ , 

as < M > c ** TZ™ ** 

9AU 1 

JJT*) 

+ 5 m c< M 


O 

■f-COS 

aAV 1 9(A° 9 A AS a A or 

«?/,* + drF** 

, dor 

4- 

9Aor 



c)(aT 

a acj 
as 






- 26 ) 


where S and n denote edge tangents and normals, A1 in parenthesis indicates an 
edge as 1 1 lusted in Figure 9 . Application of Equation (6-2fi)to the three edges of 
the triangle and equating to zero the multipliers of the terms and 

yields two equations per edge or a total of six equations for the [_c<aj unknowns 
(along edge M/ and J M , 2 = 1 , therefore the equations can be expressed 

in terms of the one variable ). Similarly, application of Equations (.e-zs) and (6 -zt) 

on the three edges and equating to zero' the multipliers of the terms and 3 

yields an additional twelve equations. The basis for these equations is summarized 
in Table 1 , With the derivatives of the area coordinates with respect to the edge 

coordinates summarized in Table 2 and the displacements defined in Equations(6~l),(c--3). 
a major algebraic effort yields the HMN equations in the form of EquationCs-i-l). These 
equations are presented in the Appendix, Section 12.5-. Application of these constraints 
yields the stiffness matrix in the form of Equation (C-^o). 
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Table 1 : HMN EQUATIONS 



57 




reprodtoibility OF ^TH t; 

OBIOINM- FMfg S 10 '®' 


TABLE 2: TRANSFORMATION EQUATIONS FOR CALCULATING HMN DERIVATIVES ON 

ELEMENT SIDES 
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The remaining degrees of freedom other than the BCIZ or HMN constrained degrees 
are then eliminated through minimum energy reduction. 


kl-Pvl-PJKJ 




where the order of [ K is to be 24 x 24, 18 x 18, or 6 x 6 as indicated 

by the options listed above. 

Since this element is based on the BCIZ element, the derivation is initially per- 
formed using the deformational degrees of freedom. The following transformation 
converts the deformational degrees of freedom into nodal degrees of freedom. 

I J “/ U, Oucj, V, @vlj t OJ 


u 2 , . , 


U-3 « * . 


J 






and the matrix [ T *J is defined in the Apendix, Section 12.6. The stiffness matrix 
is thus transformed to a 27 x 27 matrix, involving the nodal degrees of freedom. 


M-[t 
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7.0 QUADRILATERAL SHELL ELEMENT 


The quadrilateral element follows the same basic philosophy as the triangle. 

It uses the strain displacement functions of Section 3.0, and again proceeds 
from a simple basic element to obtain an element with an HMN capability. 

In this case the isoparametric element of Ahmed, et a[. Reference 3 , with 
quadratic displacement functions, was used. This element has one midside node 
per side, uses the quadratic displacement forms for the membrane and normal dis- 
placements, and introduces transverse shear deformation which are also of quadratic 
form. 

The HMN procedure basic approach is given in Section 5, The procedure differs 
somewhat in detail from that used for the triangle, due to inherent features of the 
element geometry and the nodal arrangement. 
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7.1 


ASSUMED DISPLACEMENT SHAPES 


The transverse displacement is assumed to be that of the AIZ element. Reference 

3 . The displacement assumed for the parent element has the form. Figure 10 

at 

bS = 51 F, of; 

* 7 - 0 

7 )( i * t ) 

F S =7 (i + Y)(i-?)(t-j) 

The functions F. have the characteristic that they have the value of unity at 
node i and are zero at the other nodes. 


C T-l) 


(T-'O 


The membrane displacements u and v have the same form as the transverse dis- 
placement but they include the excess degrees of freedom and nodal hotations 


a 


5 

_ V 


4 F:U: - l 6; a: - IFF F ; 0 ■ 

^ ' 


8 


V r £ /; V; ^ G.;A, £ F, 


C-T-3) 

C^-4) 
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( i-- z(t~ f)(i-r ) 

<v- (i-r)Ci-r) 

& r = -?*)(?-?') 

fi s = 7 

The functions G. are identical to the G. functions defined above except 
for G. replace t by Y) and 'O by ~t in G,. 


(7-5-j 
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The HMN functions are required to control third degree strains created by non- 
linear terms in the normal displacement. Hence, the HMN membrane displacements 
must contain terms through the fourth degree. Specifically, referring to the squate 
parent element of the isoparametric family with local coordinates ^ and y parallel 
to the sides; (1) the a displacements along lines parallel to the j axis must con- 
tain terms cubic in jf and quadratic in y ; (2) the U. displacements along lines 
parallel to the ^ axis must contain terms quadratic in * and cubic in f , which 
vanish on the sides of the element; (3) the v" displacements along lines parallel to 
the p axis must contain terms which are quartic in "j and linear in 7 . These 
three types of functions are of course repeated for the alternate functions and lines 
parallel to the 7 axis. Figure 11 shows plots of these functions. Note that type 

. r-f * v '* r 

(1) functions involve 3- a and 3- V unknowns; type (2) functions involve 2- ll 

i 

and 2-v unknowns, these being of the pillow function type; and type ( 3 ) functions 

•v ^ 

involve 4- u, and 4- V unknowns. 


The functional forms are 


Type (1) 


u — / (ny)(f-} ! ); (i-r)u-r) 


Type (2) 


Type (3) 


Ci — (t -$*■)( t-T) ; ( i-f )( y }- r l , '> 
v ; o-f) 

£->(?-!) [o?X/*- y) ( v ?) i ] ; ^ti) [o- y nHrf)C?) z '] 
<*-,) [0-7) 

n-<) [0-?)0*/)g) } ; ^h)[ 
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HMN constraints sufficient to determine these functions consist of the following: 


Constraints on 



on 

the 

lines 

*1* -1,0,* i 


on 

the 

lines 

r 

£ n 

on 

the 

lines 

3=*/ j r */ 


plus minimum potential energy relative to the type (2) functions. 
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7. 2 


ELEMENT DERIVATION 


The virtual work expressions presented in Section 4.0 are applicable to both 
the triangular and quadrilateral shell elements. The stiffness matrix was obtained 
in a general form in Equation ( 4 - 10 ^ which is repeated here for clarity 



=X /A; bo st _B x] ♦ j^g DETlllJf oil 

( 3 „ = (no,. <-ni, ;k J KH ©Jt ; , go ^ 


(7-93 


From Equation (jf'4)^ the GO matrix has the form 


l3J=tu l ,a yl ,v f ,v, 


1 < 3X l ) )&</ ) ) &*,*{ ) f , ^y, ^ 


L^J = L y : °<i ■ °< x J 



(7-12.) 

L <ji J = L u, V, cj; 6 >, s 


^3 6-4 


“3 ... 






£ * + m 

(7-13) 

U 7 , 


K J 


t J r La v a.. 


b^. b e J 


L s [a, a 3 

k 7 

Q 7 Cl £ b , b z b 3 J 

Cl- 1 ^) 


The explicit form of the GO matrix is given in the Appendix, Section 12 .7. 


67 



reproducibility of til 
ORIGINAL PAGE IS POOR 


However, since the displacement functions are prescribed for the parent element 
in the ^ ^ coordinate system, an additional Jacobian type transformation is re- 

quired to obtain the displacement gradients in the shell coordinate system,, 


L@J L , Uy , Vf 1 Vj , } ^X,« , , &y,K } 

The transormation matrix [ J ] is given in the Appendix, Section 12 o8„ 

The AO and AT matrices for the quadrilateral can be written directly from Equations 


Cl- 3.5) 
(7-lc) 


as 


if,-' (oo : - <%)&&; 




where 


L4EJ - Uf, , Afy , A K u , i ( 7 -a 

Explicit forms for AO and A1 are given in the Appendix, Section 1 1.9, 0 The curva- 

( Bendin ) 

tures (Kxx, \iyy, ftxy) are related to the bending strains by £ v 9 = £ / etc. 


8 ) 


Finally, the stress-strain relation is given in Equation ( 4 - 7 ) as 

O'; =0° -> Db u &C A 


(7-19J 


where 


L(J J = L A4 , A/ y J A/ ? /l^ , /M, , A/y, M,y J 

l a: ’J = lk\ a//, a/,;, a/; , yv,;, < , m ; , M f ;l 

and the matrix of material constants [DO] is given in the Appendix, Section 12„1Q. 


With these matrices, the gross stiffness matrix of Equation(7-9) can be generated. This 
matrix is of order 58 by 58 involving the deformational degrees of freedom Lac^J , 
the degrees of freedom to be eliminated through minimum energy reduction jLA«K fc J , 
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and the decrees of freedom to be eliminated through HMN type constraints LAo£ a 
These latter constraints have the form 


Lh](aoC a J -[ho]{a 6 j}+[hi(cj)\ {a a 7] 

L Au J = La cj, , acj 2 ,acj,, acj, , A cj r , a 44 , a cs, , Acj- t J ’■'* 

Hl;j = Hl;j k CJ K C-T-Z5) 


and the matrices [H ], [f / 0 ], and [Hi ] are developed using Table 3 . As indicated 
in Section 5 o 0, these equations can be put in the form 



( 7 - 2 - 4 -) 


and the HMN constrained stiffness matrix is then 

[k] = [c] T [k][c] It-**) 


The constraint equations are obtained by the procedure developed in Section 5.0. 
As applied to the quadrilateral, the incremental membrane strains can be written 


AC 9ACJ 3ACJ 

AE f~ if + if 


(7-aC.) 

, r _ 3AV , Scj" 3AcJ , dof 3A cjt 

^£7 ~ ^7 d*Z W + ^ 



k <3Av' 3 A AT 3*4” 

3? * c>£ 3? ^7 3$ 

3ACJ- }<jr dAcjr 

3 7 + 

(t-zs) 


69 



Table 3 : BASIS OF QUADRILATERAL ELEMENT HMN EQUATIONS 









Ten HMN type equations can be written by evaluating these strains along various 
edges or locations within the quadrilateral and equating to zero the multipliers of 
terms in the strains,, The basis for these equations is summarized in Table 3 „ 

With the displacement functions defined in Section 7„1, an algebraic expansion of 
Equations ( 7 - zc - 2 . 8) , yields the HMN constraint equations in the form of Equation (? -'2.1 
These equations are presented in the Appendix, Section l2 0 N. Application of these 
constraints yields the stiffness matrix in the form of Equation ( 7 -Z50. 

The remaining degrees of freedom, [<X^ ], are then eliminated through minimum 
energy reduction, 

where the order of [ K*” ] is 40 by 40, five degrees of freedom at each of the eight 
nodes. 
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8.0 COORDINATE UPDATING TRANSFORMATION 

A special type of updating is used in the nonlinear solution procedure. This 
updating essentially uses convected, or material, coordinates, pertaining to 
the element start-of-step orientation, to set up an updated Lagrangian form- 
ulation of the incremental equations. The updating requires a special trans- 
formation to permit the entire solution process to be purely Lagrangian in 
character. For each new incremental solution the accumulated total and the 
incremental forces and displacements are referred to the end-of-step element 
base plane coordinate system of the previous step. The updated end-of-step 
element coordinate systems are determined for the iteratively converged shape 
of the structure at the end of the particular increment in question. Strains 
and stresses are computed with reference to the updated systems permitting a 
determination of total strain and stress values without the necessity of summing 
increments. Updated coordinate systems are also used as temporary coordinate 
systems to facilitate the iterative calculations. In this case the 
purpose is to obtain an exact strain/stress state on which to base the calcula- 
tion of residual loads and the iterative equilibrium corrections. The sequence 
of computations is sketched briefly below: 

o Determine end of step element baseplane coordinate system (converged), 
o Compute stiffness matrix, loads, etc., referred to this as a new 
start-of-step system. 

o Compute incremental displacements, referred to above, 
o Determine temporary element end-of-step coordinate system, cor- 
responding to computed incremental displacements, 
o Transform displacements to the temporary coordinate system* 
o Calculate strain and stress (total referred to the (temporary) 
end-of-step coordinate system**) 
o Determine residual loads and check convergence 
o Correct incremental displacements, to eliminate residual 
loads (uses start-of-step stiffness matrix) 


* This is a special transformation, to be derived in this section. 

This is not by summing increments; It is a direct total calculation. 
The strains and stresses are of the Lagrangian type. 
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o Combine correction displacements with prior displacement 
increment and repeat determination of temporary end-of-step 
coordinate system, strain, stress, and residual load cal- 
culations, etc, 

o Identify the last temporary end-of-step coordinate system an 
the final, converged end-of-step element coordinate system 

At this point the entire calculation passes on to the next loading increment. 

The remaining small residual loads are carried on into the next step to avoid 
accumulating error. A more detailed description of the solution procedure 
is given in the next section. The above brief listing serves to explain the 
role of the special displacement updating transformation which is employed to 
compute total strains. The derivation of this transformation is the primary 
purpose of this section. 

Consider Figure 12. An element is shown referred to three rectangular 

cartesian coordinate systems: the initial one x 1 for the undeformed structure; 

, -i H -2 ° 

a coordinate system x whose x , x coordinate plane is parallel to the element 

fV»j Oj] ^2 

in its deformed state; and a coordinate system x whose x , x plane is parallel 

to the el ement -after a further increment of displacement has occurred. For 

simplicity the figure is drawn in only two dimensions. The derivation given, 

however, includes all three components of displacements and coordinates. The 

original coordinates x' are a cartesian version of the material coordinates. 

o 

In the cartesian sense, x 1 can be considered to be convected with the deformation. 
_i m ° | 

Thus each of x and x have Companion, col inear x q systems. As the deformation 

progresses, the material coordinates are deformed and rotated but still partially 

coincide in direction with the x and x coordinates, as follows: the convected 

x^ axis coincides in direction with the x' and axes and the convected x^ - x^ 

° — 1 — 2 'v 1 ^2 o o 

plane coincides with the x , x and x , x planes. In addition, the side of the 

element containing nodes 1 and 2 coincides in direction with the convected x , 

— ] , ° 

and x , and x axes, and node 1 is at the origin of all coordinate systems. 

— ] j 

Thus, the coordinate systems x and x are the rectangular cartesian counter- 

parts of the element convected coordinate system, with the conventions of 

coincidence of the origin, the number one coordinate line, and the 1-2 plane. 

Note that the metric of x and x contains only zero and unity, while the metric 

of x 1 does not have unit values after deformation has occurred, 
o 
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Displacements, both total and incremental, used in each step, are referred 
to the current coordinate systems. Hence for the increment whose start-of- 
step system is x 1 , the displacements u" 1 (x^) , AH' (x^) , and u" 1 (x^) + Au" 1 (x^) 
are components in the directions of the x 1 axes. The total strains and the 
incremental strains of the Lagrangian type, referred to the x 1 system, viewed 
as the convected x^ system, are derivable from these displacements provided 
that they represent, in the appropriate coordinate system, precisely the 
deflections from an undistorted shape. At the end of the step, the u 1 + au' 

f\, j 

are transformed to refer them to the end-of-step x system, again as deflec- 
tions from an undistorted shape. The derivation of a transformation of 
IT' + AlT 1 from the x" 1 system to the end-of-step x‘ system such that the re- 
suiting displacements u adhere to this definition is given below. 


The undeformed shape in the x^ coordinate system is described by the position 
vector (.Figure 12) 



• 


x ! G. 
o i 


(8-l) 


where here G. are the base (unit) vectors of the x q system. Here x q are the 
coordinates of points on the undeformed shape of the element. If the unde- 
formed element is placed into the x 1 system (shown dotted in the figure), 
its position vector relative to the origin of the x' system has the same 
form , 



i 

x 

o 
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Figure 12. ; cooroihate -systems ano befikitions 
coordinate system updating procedure 


For 
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where here g. are the (unit) base vectors of the x 1 system. The actual shape 
of the deformed element is given by 


P = (x + 
o 


-* 

9 ; 


(*-sJ 


where here the u 1 are implicitly defined by this equation to be the deviation 
of the deformed element shape from the undeformed shape, referred to the x 1 
system. After the incremental displacement, 


i , —i 


P + AP = (x - + 
o 


+ Au 1 ) 


• 4 - 

9 ; 


( 8 - 4 -) 


This vector, referred to the x 1 system, actually defines points on the element 
depicted by the solid lines in the figure and "attached 11 to the x' system. 

That is, the Au 1 are the displacements which carry the deformed element of 
the x 1 system to the deformed element of the x 1 system; x 1 are the end-of- 
step coordinates of the step for which x' are the start of step coordinates. 

The position vector of a generic point of the element in the x' system is given 
by 

? = c < * cr ) c?- 5 ) 


where IT. are the (unit) base vectors of the x 1 system. The implicit definition 

r\j j 

of u as the deviation from an undeformed shape is again apparent. Note that 
u' defines the deformation condition which corresponds to u' + Au'. 

— j (\, j — — 

Let R be the displacement of the origin from x to x . Then P + AP « R + P, i»e M 




■>n = 


R + (* 0 + “ J b. 


(8-c) 
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Def in 


l"g V> ; • >. = X; 3 


( 8-71 


= R. • V\ - + * Q 4- 


This gives the desired transformation and its derivative with respect to J 

the material coordinates. 

!= ’ 

( 8 - 5 ) 

i- 

(S-lq) 

In the above, distinction between contravar i ant and covariant components 

has not been made, since the coordinate systems are cartesian. 6 

These equations transform the end of step deformations referred to the x 1 

f\, f 

system to the corresponding start of step deformation referred to the x system. 

This transformation is used for updating the derivative freedoms of the 

‘ * 

element in order to define the total deformation state of the elements both 
at the end of the step and equivalently, for the start of the next step. 

The displacement freedoms can be transformed using the above equations, or, 
alternatively, they can be. obtained by simply subtracting the deformed nodal 
coordinates from the undeformed ones, in the x' system, at the end of the 
step. Note that in the case of the triangular element the above are "total" 
derivative values rather than values of the BC1Z type which omit the simplex 
type terms. Thus, the BCIZ transformation is required to be applied to the ? 

above types of derivatives in order to compute strains in the element using 
the standard BCIZ formulas. 
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9.0 SOLUTION PROCEDURE 


Nonlinear solution procedures and the formulations of nonlinear finite 
elements are interrelated. Frequently derivation details are specifically 
tailored to planned solution procedures. In addition, precisely how the 
element equations are used, including updating, transformation, summing 
increments, etc., in the solution procedure has important implications re- 
garding the performance of the element and the validity of the element 
derivation steps. The elements developed in this research haveclose inter- 
relationships between derivation details and solution procedure steps. This 
section discusses the solution procedure and related element derivation 
deta i Ts . 

The solution procedure is a modified incremental Lagrangian type. Each in- 
cremental step is subjected to iterative equilibrium corrections which in- 
volve the calculation and reversed application of residual loads. To accom- 
plish this, element total strains are computed directly (not by summing in- 
crements), corresponding to the total displacement state. The total strain 
calculation employs an element coordinate system updating which involves the 
special transformation described in Section 8.0 The overall procedure can 
be classified as either a total or an incremental Lagrangian type, with up- 
dating . 

The steps in the solution procedure are described below: 

1: The stiffness matrices are generated based on the previous end-of-step 

displacements, stress state, and updated element coordinate plane. The 
element stiffness matrices are transformed to the solution coordinate 
systems, using the locations of these systems for the end of the pre- 
vious step. The element matrices are then merged to form the overall 
stiffness matrix, referred to the solution coordinate systems. (This 
is not the global system.) In iterations pertaining to a given load 
increment, the stiffness matrix and related transformation matrices 
a re not changed . 
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2 . 


3. 


k. 


The load increment, including the remaining residual load from the 
end of the previous step, is obtained and transformed to current 
solution coordinate systems. 

The incremental displacement is calculated. This is the first 

estimate, uniterated, of the displacement increment for the current 
load increment. 

The above are transformed into the element baseplane coordin- 

ate system 





5 . For the triangle only, for use only in a later calculation of the 
minimum potential energy freedom incitements, the above are 

transformed to the "deformational" incremental displacements. These 
are the increments used in the BC1Z element derivation. 

For the quadrilateral this would be an identity transformation. 


6 . 


The element incremental degrees of freedom are summed into the pre 
vious totals to form new total element freedom values. This uses 
the values ' from the previous incremental calculation, as 

defined in Steps \k and 15* 


* 0 > 


*0-0 Ci) 


The sum obtained in this way contains the current estimate of the 
displacements of the element, relative to its undeformed shape, and 
referred to the current start-of-step element baseplane coordinate 
systems. The displacement derivative freedoms obtained are total 
rather than deformat iona 1 , That is, for example, • includes not 
only the BCIZ deformational contribution, but also the contribution 
due to the nodal deflections. Since the latter are the "simplex" 
values, we have used the subscript "SIM" to denote the fact that the 


t 
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gSSKSSM* 


7. 


8. 


derivative freedoms_ in S contain the simplex contributions. 

5 lh^ 

The incremental deflections are transformed to the global coordinate 
system and summed 

C3> 


[tr M \ - fu a "^ «- l ^ Ci> ] 


Here the are the nodal deflections contained in the displacement 

^ from Step 3- 


The element geometry and location in the global coordinate system is 
updated in order to update the element baseplane coord inate system. 

- {x . ci -° ) * 

Note that this updating is not a generic element updating. See Sec- 
tion 8.0 


9. A new matrix of direction cosines relating the updated element base- 
plane coordinate system to the global system is calculated. These are 
called 

IT (J) ] 

[T^] is defined such that nodes 1 and 2 are on the updated x axis 
and either (for the triangle) nodes 1, 2 and 3 are In the xy plane, or 
(for the quadrilateral), the nodes are displaced from the xy plane 
such that their RMS departure is minimized. 

10. The current nodal coordinates referred to the updated baseplane co- 
ordinate system are calculated. For the triangle, the set of values 
includes only x^» x^ and y^> due to the definition of IT^] given 
above. For the quadrilateral it includes all x. , y. except y^, yg, plus 
the z coordinate for all nodes. The equation is 

[/*>]■» [t‘ w 1\.5: ch 1 

Here, .IX contains the quantities X^f’n where X^^' r. = X^^: - X^' 
refers to the X, Y, and Z global coordinates. 
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The result of this calculation is that the nodal coordinates are known 
for the current state of deformation Tn the current end-of-step updated 
element baseplane coordinate system. These data permit calculation of 
nodal deflections relative to the undeformed element referred to this 
same coordinate system. The updated element coordinate system is used 
essentially as a material, or convected, coordinate system. 


1 1 . 


The element total nodal translations are next calculated from the 
equat i on 


C ~ C 3)7 


c 15) 7 

1 * 5 

— * 

r \ 



The notation ( ) denotes the fact that these translations are referred 
to the current end-of-step element baseplane coordinate system. The 
u VJ include the in-plane defections only. For the triangle, the normal- 
to-the-updated-baseplane deflections are equal to zero. For the quad- 
dri lateral, the normal deflection are set equal to the changes in the 
normal coordinates determined by the RMS fit of the four nodes to the 
updated element baseplane. 


12. The incremental minimum energy freedoms are computed by back substitu- 
tion using the same constraint matrices used originally to constrain 
the stiffness matrix referred to in Step 1. This calculation effective- 
ly uses the tangent stiffness of the element for the current increment 
referred to its s tart-of-step coordinate plane. The computed quantities 
are 


l . 


<sv 


„ ti) - 


| 

i 1 . 

*<*. 

•»* 




depending on the particular HMN constraint options employed. 
13* The minimum energy freedoms are next summed to obtain totals 


R'l 


oL, 


Ci- 0 




G) 


These quantities are not subjected to any geometrical transformation 
type of updating. 

]k. The displacement derivative freedoms are next transformed to the end- 

4 

of-step element baseplane coordinate system using the R transformation 
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of Section 8. The equation is 

f 


~*G> 

o 


Sim 


S ■ * f«*5 


where the notation ( ) denotes displacements and differentiations 
referred to the end-of-step element baseplane coordinate system and 
the subscript "SIM" again notes that the derivatives obtained are 
total rather than BC!Z deformat ional quantities; they contain the 
simplex contr ibut ions . The transformation is applied selectively to 
the derivative freedoms » since the nodal translations 

themselves were computed in Step 11. 

15- The nodal translations from Step 11, which are also referred to the 

end-of-step system, are entered into the46 c ,y f vector, replacing 

^ ^ ^ * ( i ) ? 

the corresponding translations term by term. At this po i nt J 

describes completely the element deformation state, relative to the 

undeformed element shape, and referred to the end-of-step element 

c** * ( j ) -, 

baseplane coordinate system. Note that the resu 1 1 ing | S s j f is the 

vector to be used in Step 6 in performing the next incre- 

1 SIM 

mental calculation. £$ s | M j ' s 5avec ^ f° r This purpose. 

16. The BCIZ transformation jj*3 is applied toS*.^ ' to remove the 
simplex contributions to the derivatives. 

[ s‘™} = 

r ] \ -j 

This prepares the vectoriS J 1 for use in total strain calculation 
in the end-of-step coordinate system. The resulting strains will be 
total Lagrangian values. 

The total HMN f reedomsjot^ ^ are next calculated. These freedoms are 
determined by constraining the polynomial forms of the strains along 
the element sides. The equations governing the total HMN freedoms, 
in the present form of the computer program, are of the form 


17. 


v D J 


C 6 *^ 

A. 

y ^ 






for the triangle, and 
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* iV£*r>] P"} 


for the quadrilateral. Here W are the HMN function amplitudes and 
are t * ie m ' n ' muni P° tent ' a l energy parameters .[ 0 and 


© " (j) ^ are nodal values of the derivatives of £ w ^ J and £w ^ ^ J , 
respectively, with respect to the current x and y element baseplane 
coord i nates . ^ ^and|© w are contained in^S and 

^ 8 (^)j where these vectors are as defined In Steps 14 and 1 5 • {_ S ^ 
refers here to the deformations at the end of the (j-l)st step, re- 
ferred to the end-'-of-step baseplane of that step. j^S J ( refers here 
to the deformations at the end of the jth step, referred to the end 
of step baseplane of that step. Equations similar to the above are 
used in Step 1 to constrain the stiffness matrix for the incremental 
calculation. From these computations the matrices [60^] and {6 1 ^ ^ 
are saved. These saved matrices are used to compute the total HMN 
freedoms at the end of the jth increment. 

The use of£# 1 ^ j and£©^ together is an approximation which has 

been found of small effect. It can be made as exact as desired by the 


use of zero load increments. The exact equation, of course, would 

result from reforming the matrix using the currents 

mnp w w 

referred to the end-of-step baseplane coordinate system. 




For the triangle, 
o) 


{< l 


u) 


ti) 


\ l V,uU ^ 


according to which HMN option is chosen. 
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!t Is noted that the above HMN procedure, particularly if employed 
with occasional zero load increments, determines the HMN freedoms to 
control the total Lagrangian strain, obtained by reference to a cur- 
rent, updated, element baseplane. Thus, the HMN freedoms are not 
subject to accumulated error. 


18. 


At this point the current totals of all freedoms are available, re- 
ferred to the end-of-step element coordinate system. These totals 
consist of two sets, as follows: 



The former takes^fi $ from Step 14 and is used for summing to obtain 
totals for subsequent increments, as described in Step 6 and further 
explained in Steps 14 and 15- The latter takes^S from Step 16 

and is used in the jth increment for HMN function calculation, just 
described, and also f or strain and stress calculation. The latter 
calculations permit the determination of residual (error) loadings, 
and hence, iteration of the solution to improve overall equilibrium. 


Steps 19 through 22 are calculations performed for each integration point 
of the element, leading to a virtual work integration performed in Step 
23- 


19. The matrix of displacement gradients is calculated at all integration 
points 


& 



Go 


v», v\ 


9L 
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20. The strains are calculated 
_ Ci) 


= C A °. 


+ — AG *0^^^ 

t. > 


v\ 


These are the true Lagrangian strains, referred to original material 
"fibers" provided the membrane deformations are small. 


21. The stresses are calculated 


cr 

W\ 


_ cr 


DO 


CiT 


22, 


Quantities are next calculated which will provide virtual strains at 
each integration point and lead to the virtual work integration 


23- 


** G) 





(.3) ^ O ) 


~ cr a, 

w v\ yn 


~ (wt) QO 


M 


^ * G) 


^ G) 

pVA ^ p 

Here (WT) are the weighting numbers for the particular integration 
scheme employed. The stresses are taken from Step 21 and describe 

P i) 

the current end-of-step state of the element. The subscript in Q^ J 
denotes the particular element freedom corresponding to which the 
virtual work of the 0^^ is sought. The result of summing the 
over all integration points is the total internal virtual work cor- 
responding to a unit virtual increment in the m^ freedom, and pro- 
vides the current generalized load on the element corresponding to 
that freedom, referred to the end-of-step coordinate system. 

The current generalized loads on the element are computed by summing 
over the integration points 


- (.*>) 
Q 


RNTtpN 

Points 


24. The constraints associated with the HMN and minimum potential energy 
conditions and, for the triangle, the BCIZ transformation, must be im- 
posed on the Q ^ in order to obtain generalized forces which correspond 
to applied loadings in meaning and ordering. The equations are, respect- 
ively 


85 



REPRODUCIBILITY OF THU 
ORIGINAL PAGE IS POOR 


i - i* p l «“'] 

r ifcG)? <■ . C 5 ) 7 t* "1 r i> t3 ^ 1~^ C r 

[a J = ivS*I K 5‘iH K w] l 01 1 J 

[S c5> ^ = [T* <i> l T [ ta* t5) } 

l 

Definitions of terms in the above equations are given in earlier, 
related sections of this document. The j are the element nodal 

forces referred to the end-of-step element coordinate system. 


25. 


The element forces are next transformed into the current, end-of-step 
solution coordinate system, using the updated matrix IA^ ] (see Step A). 


1^} = lX U) ][S a) } 

26. The \Q. }for all elements are merged into a single vector, each member of 
which refers to the appropriate end-of-step, nodal, solution coord- 
inate systems, and the residual loads are computed by subraction from 
the applied loads, also transformed to these coordinate systems, 




- f £ c5) l 

l_ fvmisp} 



27- A convergence check is made to ascertain if the residual loads are suit- 
ably small, If so, the solution procedure returns to Step 1; If not, 
the steps below are carried out. 

The Iterative calculations described below are repetitions of the steps 
described previously, using different starting values of as 

described below. In the current form of the computer program, the stiff- 
ness matrix for a given load increment is not updated for the iteration 
of that load increment. Of course, updating can be accomplished simply 
by the use of a zero load increment. It is again noted that the residuals 
remaining after convergence has been achieved for a given increment are 
always added into the loads for the next load increment. 
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28. The incremental load is set equal to the residual calculated in Step 26 . 

2 9- Step 3 is repeated, using stiffness matrix from Step 1. The results 
are correction values of£&q^j. These are added directly to the 
previous results of Step 3, the result becoming an improved estimate 
of ^^q . Various schemes to accelerate convergence have been us^' 

here, mostly based on constraining the bending type corrective incre- 
mental displacements during the iteration. The constraints are applied 
on every other iterative step. The best procedure found is to use the 
membrane components of the stiffness matrix from Step 1 in every other 
iteration, and the full stiffness matrix in the alternate iterations. 
This cannot be recommended as a suitable procedure for general usage, 
however, as its success may be specific to the particular problems 
under study. 


At this point Steps 4 through 27 are repeated exactly as described 
previously. The updating now is to an improved end-of-step element 
baseplane, and all transformations refer to this improved coordinate 

^ * C j - 1 ) 

plane. In Step 17 the HMN calculation again makes use of the6? w 
the values from the previous load step, not those from the previous 
iteration of the current load step. 



10.0 RESULTS AND CONCLUSIONS 


Pilot computer programs for both the quadrilateral and triahgular shell 
elements have been written. These programs can accommodate small assem- 
blages of elements, and each includes a nonlinear solution procedure 
adapted to the element formulation. 

Checkout of the computer programs has included successful compilation, a 
limited unsuccessful numerical application for the quadrilateral, and a 
detailed investigation of two different problem solutions for the triangle. 
The triangle has been used to work out solution procedure difficulties and 
to investigate the effect of the HMN formulation in nonlinear problems. 

This work has been successful, and it is felt that the value of the HMN 
application to two dimensional problems has been demonstrated. The fol- 
lowing discussions of results pertain to the triangular element studies. 

The problem solutions obtained serve to check the membrane only and the 
coupled membrane-bending behavior of the the triangular element. The 
various transformations, updating, etc., in the solution procedure are 
also checked by these problems. The problem solutions themselves are dis- 
cussed in Section 10. 1 and 10.2. The difficulties involved in the solution 
procedure are discussed below. 

The solution procedure includes the calculation of residual loads and itera- 
tion to improve the structural equilibrium state. The basic solution pro- 
ceeds on a stepwise linear basis. In general, the stepwise linearity gives 
rise to large residual forces due to the effects of the rotations. It was 
found that divergence of the solution occurred when the residuals exceeded 
certain magnitudes, in particular, magnitudes approaching structural In- 
stability loads,, These limiting magnitudes were obtained with quite small 
load steps, with the consequence that practical solutions of large deflec- 
tion problems required an excessive number of load steps. This difficulty 
has been resolved, and the solution procedure now performs well for load 
steps of essentially arbitrarily large size. Convergence is extremely rapid 
and can in most cases be obtained in a single corrective iteration. The 
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current solution procedure Uses a constrained, membrane only, stiffness 
matrix for alternate iterative corrections, and the total stiffness matrix 
for the complementary alternate iterations. 

The quadrilateral and triangular elements both permit representation of 
either flat plate or curved shell structures. Curved shell checkout pro- 
blems have not been attempted. It is noted, however, that the equations 
by which the elements represent nonlinearity due to developing rotations 
and curvatures for flat plate problems are the same as those which account 
for the effects of initial curvature. Thus, the flat plate checkout work 
has effectively dealt with the curved shell features to a significant extent. 

10.1 TENT PROBLEM 

The tent problem was designed to provide a simple check on the program 
flow and solution procedure. The loads and geometry are such that no bend- 
ing occurs. The flat membrane stiffness, stresses, and most basic trans- 
formations and matrix operations are checked, however. The basic configura- 
tion, idealization and boundary conditions are shown in Figure 13. To re- 
duce the problem size a plane of symmetry was assumed parallel to the Z axis 
and passing through nodes 2 and A. This problem was successfully solved 
as shown by the results given in Figure lA. The seven step triangular 
element nonlinear solution is exactly the same as the theoretical solution. 

10.2 BEAM PROBLEM 

The beam problem is a simple example intended to exercise both bending and 
membrane behavior and to check the performance of the HMN procedure. The 
structure is a three element, initially flat plate shown in Figure 15 - The 
figure also shows the local coordinate systems of the elements and the struc- 
tural constraints imposed. As described in earlier sections, several options 
are programmed for the HMN procedure. The one used in this sample problem 
solution constrains the high degree membrane strains parallel to each 
element side, along the particular side in question. There are thus six 
HMN unknowns per element. The HMN functions provide displacements tangential 
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to a side, of quintic and quartic form in the side-length variable. All 
HMN displacements normal to the sides have been set at zero. All HMN 
interior functions have been retained and controlled by minimizing the 
total potential energy. This HMN option is option #1 as described in 
Section 6.2. 

The overall stiffness in bending of the cantilevered beam of Figure 15 under 
end load is about 15% too high for the initial linear behavior. This is 
caused by the coarseness of the finite element modeling and the displacement 
function of the basic BCIZ elements, specifically the constrained BC I Z 
cubic displacement form. Figure <16 shows the end displacement versus load 
results, which are almost linear through the load range considered. The 
deviation from linearity seen in the figure is consistent with the effect 
of transforming the applied vertical load to the element solution coordinate 
system. The insert in the figure shows the distribution of bending moment 
on the edge of element #1 over the length of tne beam. This distribution 
supports the conclusion that the element arrangement is too coarse for 
application to this problem with the basic BCIZ element. The total nodal 
moment reaction at the support end of the plate is in equilibrium with the 
app ) i ed load . 

Figure 17 shows the spanwise stresses along side 1-3 of element number 1. 

The successive plots show the stresses at three load levels and both before 
and after a zero load step. it Is recalled that the zero load step is used 
to allow the HMN procedure to improve In accuracy by implementing the most 
recently updated nonlinear strain values. For the lowest load level shown, 
the diagram on the left corresponds to zero values of the HMN function; no 
HMN participation has yet developed. The diagram on the right corresponds 
to essentially complete HMN effectiveness. Each of the two diagrams on the 
left for the higher Ipad levels result from partial effectiveness of the 
HMN functions, while those on the right again correspond to essentially 
complete HMN effectiveness. It is seen in each case that, the HMN procedure 
has reduced the stresses significantly. In addition, the forms of the 
residual stresses after updating are of relatively low order, apparently 
quadratic, of the BCIZ type, rather than of the high order 
form which would result if the HMN procedure were not used., *-;» .onch 
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case the higher degree forms, which are seen to he predominantly cubic, 
occur only prior to the zero load step, when the HMN functions are not 
fully operative. Similar calculations have been made for sides 1-3 of 
elements // 2 and //3. The stresses parallel to the side were computed 
before and after the zesvj load step. The results were essentially the 
same as those shown in Figure 17, with again the HMN functions reducing 
the stresses to very low levels. In this case the remaining residual 
stresses were distributed linearly over the length of the side. As 
would be expected, the stress parallel to the side has identical values 
in the adjoining elements. Finally, all components of the stresses in 
the interior of the element are very small also, of the same order as those 
parallel to the sides. 

As another check on the levels of the membrane stresses, the membrane 
strain energy was computed and its magnitude compared with that of the 
bending strain energy. This check showed the membrane strain energy to 
be negligibly small compared to the bending energy, and further verifies 
that the element’s control of the nonlinear strains is functioning well. 

The study of the stresses along the sides of the element has shown 
(Figure 17) that the residual membrane stresses are of forms which are 
necessarily produced by the BCIZ membrane functions and the lower degree 
nonlinear strains due to bending. The higher degree nonlinear strains 
due to bending have been eliminated by the HMN functions. This raises the 
somewhat academic question of why the BCIZ membrane functions have not 
completely eliminated also the lower degree nonlinear strains due to bending 
since it would appear that the lowest energy state would have zero values 
of these strains. A complete answer to this question is not known at the 
present time. It appears most likely that the cause lies in the use of 
HMN option //I ; HMN option #2 would provide a more complete control of the 
stress state. 
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12,0 APPENDIX 


1 2.1 TRANSFORMATION MATRICES 

1 2. 1 . 1 Triangular Element 


Global to Local Transformation 


Input: 


Define 


Nodal Coordinates 
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Then the coordinates of the nodes in the local system are 
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where the matrix of direction cosines are given below. For simplicity omit the 
superscript } denoting the step, i,e., let 
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Then at any step }, the matrix of direction cosines is recomputed. 
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Global to Solution Coordinate System Transformation 


At a node, average the unit normal vectors k of all M elements having 
this common node. Figure 2 . Again it is implied that the calculation will be re- 
peated for each solution step 
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Since i and j are arbitrary, they are selected to be in the XZ and YZ 
planes respectively. 
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Solution fro Local Coordinate System Transformation 


— .... — ,i ■ .■ ' — . .. ■ 

From Equations (l'Z-.A) and Q/z. i?) f for each node of each element at a given load increment 



The displacement gradients transform as a second order mixed tensor 
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Expanding this equation and using Equation the transformation matrix for dis 

ment gradients for each node of each element at a given load increment is 
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where 



A 

The derivatives with respect to - X 5 have been deleted from the stiffness matrix. 
The matrices given in Equations l£1.ili)and Q.2.18) can then be re-ordered and assembled 
for each node to give the transformation matrix 
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12.1.2 Quadrilateral Element 

Mean Plane Calculation and Global to Local Transformation 
Input: Nodal Coordinates 

oc\r\zr) 


Node 1 



A mean plane can be shown to exist that passes through the midpoints of the 
straight lines joining the corner nodes. Figure 1 . For simplicity, the superscript 
denoting the step is omitted since the mean plane and the matrix of direction co- 
sines is recomputed at each step. 
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The matrix of direction cosines is 
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(1*2 .A 


(12- . 2 l) 


(12.22) 


(.12. 22.) 
(tOUT.) 

Cl 2. 23) 


Cl 2 .24) 



(M&m&h PAGE IE PGwi. 

The coordinates of the nodes in the local coordinate system can then be calculated 
using the matrix of direction cosines and the defined origin of the local coordinates 

M (x-x.1 

U;) (%;'Zo) 


where the Z coordinate is the initial shape w . 


Solution to Local Coordinate System Transformation 


The translational degrees of freedom transform like the triangular element, Equation 12.13* 
The rotational degrees of freedom transform as vectors, so that at each node of each 
element of a given load increment the transformation has the form 


/J US 


where rotations about the 


matrix. 


Z axis ( A©,) ^ ave been deleted from the stiffness 
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1 2.2 


TRIANGULAR ELEMENT GO MATRIX 


The area coordinates are related to rectangular cartesian coordinates (Figure 9 ) by 



<jz /ZR 

e /E ' -y,/2^ 

e„ - (x,-xJ/2# 

& 22 " f/j X? ^ 

- -x 3 

- K/zff 

£ zl - £j t - ^ = O 
Then the GO matrix is given by 


GO 


^ /£ 

GO^ 

- 



= 


G.O rii 

r 

if. 

dX 

GO s-,9 

- 

df*. 

TT7 

GO XT, ,7 

- 

IT* 

GO SilS 

- 

7T 

105 




(li.lB ) 


(12 .29) 


(l2.1.o3 



REPRODUCIBILITY OF THB 


v TMfti] Tg p Q , 


ao. 




^26 " JX 


y, Z7 

= 

ax 


*- /a 

^6,/c 

^23 


= <?*, 

&°6 t 8 

- Jfj 


dF* 

dy 


° ; f 


e)Fy 

v ^T 

^^6,26 

a-e r 


- -3T 



d f F, 

'7/* t 

ao 7i , = 

Ji'ii 

a** 

Go 7 „ 7 - 

av? 

c)x a 

G 0 = 

c3 ^ Pv 

ax*- 

GO, „ - 

aw 

8x z 

G0 7j n -- 

<^ Z F 0 
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ao, jB 

Ff, 

z 

G.Og^ 

f££L 

A A 

<^Fj 

^8, i? 

~ cty 1 

GtOg.i? 

- av* 

Jf 

6o s , It 

d*-F* 

- 

&Og,vi 




«?v, 

z <Wy 

S0, A 

<5 l X t 

= “373^ 

ao V7 

_ 


dO^ )( g 


&\tc 

cJ'Fs- 


&0, |8 . 7 

<**Fc 

C?* 0>JJ 


o" 

~ fi0r ( T 


'- fio, g 

fiOi,a 

- ^ /0 
~ ^OV.ig 

d.0 /,<» 

0>4-, 

o>* 

r,Z£ 

„ <3<*v 

c>X 

dO JjVC 

^ F,< 

- Cos*, <* x 

i,m 

^ 4 T 

= £<**« ^ 


fi0„„ = 

^5,tC 


fio„, -' 



&0 (lt .= 

60** 





^°/,3= = 

PGtf" 

d* 




C^Mtl 


CoSoC , 

3x 

£<W 

dot 0t 3 

«r^ ( 3 


do,,* - 

&0 i,u " &^sr,n 


Q> 

r 

ll 

&0r ( t7 

0><*5 

G0 h <H 

dx 

Go h ^~ 

C?<*0 

5T 

Go,^t “ 

Mt 


6^.n " 

c/Ho 

e°i *3 37“ 
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GO i t 3H ~ 

m sm <*, ~jy~ 

Go^i5- - 

- j-/v? 

Mti 

TT 

&u= 

^ Ms 

- 5in&Lz_ c ,)% 

GO , i37 = 

-S/no( z 

G 0 , ,xs - 

- 5//? &c 3 

"$ rt & 

tst 

G<9 / 3 i" 

, 3tfa* 

-Stn<X 3 

So 2< , r 

GOt.j 

6 ^ 2 , , 0 = 

GO^tc, 


* 

G 0 «.. r 

QOzz. - 

GO 6,8 

G 6 2rJ - 

G0c,9 


go 2> „ - 

GO<,,n 

GxOz'H — 

G>O w% 

G 0 2|W - 

GO 6,26 


Go iZi s= 

GO c,,z7 

GQ.vo = 

G&Z'ZQ ^ 

3& i 
dy 

aov 

' GO gHZ - 
G 0 Z 3O - 

SGz 

~dj~ 

9 Gf 


G 0 ?vy - 

GOzjz ~ 

d&i 

96 0 

J(j 

GO^-- 

<9 Hit 

tosoc, 

Go Z c tl - 


Go z ^ 

CosoCj, '9 lj 

GO ZiH t ~ 

Q, OS 0 (.£ Jy 

Go Z)SO - 

C-OS^i^. 

GO «* 

J//u 

c °sd 3 

G 0 Zl 3 y- 

Go Zl » -- 

3 H n 

-Smct , 

9 Hz.z 

-StnoLz 5y 

Go zar ~ 

- Ga^ 

-sin °6 
~-S(n o ( 3 

9 Hzt 

■3T 

3 // 13 

G 0 ? , y6 = 
GG*,« 55 

. , 

"5m o <: 2 - 35 - 
1 / AMh~ 

~Sih. oCjj 

GO,,. = 

GApr 

GO 3, i3 = 

GO r //4 


G<V * 

GO ^r,£a- 

GOi,r - 

GOV,* 

GO** - 

G o r# i 


= 

60^77 

GOj,»- = 

G <9^-,/# 

GO 3,2, - 

GO^u 


G0 3/ti * 

60^27 

GO j,*/ — 

60,,y o 

Go^v, — 

GO /jHZ 


& 0 3 pr~ 

GO/,yy 

GO ' 3|Z < = 

<30,,*. 

GO 3j3 / “ 

GO t t 3e> 


GO,^ = 

GO /,32 

GO J)V .~ 

G 6*7,3 y 

GO 3 t *t 7 ~ 

-Go i.str 


G Oj vg - 

~Go ,,3<* 

GOj,^= 

go /;37 

GO 3 / y 0 = 

-G 0 Ii39 


G 0 3f , - 

-GO/^t 
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* 

M 


G0 33r = 

GC> //vo . 

- 

(SO,,« 

G0 3 l 3*f - 

GG/,y r 

GG» 3(M = 


GO 1 . 1 , - 

G O ,'f, 

— 

gg., 7 

ea,* - 

GO,.jc 

GOij'Zi - 

GO^, Z r 

GOi, r = 

G0 6j « 

G0y )6 = 

60vi 

GO^,n - 

GO (,,a 

G^V = 

GO Cilt 

fi^u s 

G,Oc,zc. 

GO %z i - 

GO 6 |2 .7 

G0i tV/ - 

G0 Z ^ 

GO%11 ~ 

G 02, y2 

GOiH<~ 

G 0 g^yy 

G£y,zi - 

&0 Zi zt 

G.0^'3, - 

G02,jo 

G0<4,JJ ~ 

GOz^z 

~ 

~G0 Z3 ^ 

G 0*1,1 7 “ 

-G0«^ 

&0 % * = 

- Gi9 2/ 3* 

Go^ii ~ 

~GO z> 37 

GO** = 

-GOz'jt 

G<2y, r , - 

"GCg, 3 <? 

Go^ - 

Go,,, 

G-0*t,3X ~ 

GO e> n 

G Gy,jr - 

G0z,v« 

G J7 — 

G0g,lj 

GO«,n = 

G0 Z ^ 

G0 1,37 - 

GOz t n 


*v 
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TRIANGULAR ELEMENT AO & A1 MATRICES 


1 2,3 


ftOsr - 
/?O ai = 
RO&(, ~ 

^ s 


1 

ftliss ' 1 



■^T 


1 

“ -2 

9as° 


-7T 

i 

/«» = 1 


51 1 



- 

1 




ft < 03$- 


« 0 

0OJ 

IT 




hd 3( . 

■=* 

9ar° 

7T 




no* 

s ' 

Z 

7** 

ftl 

//J 7/7 

= -i 

/7/ V7/ - 

POk 

« 

#v° 
“ ^ 

«** 

= -; 

>W **3 = 

ft 

- 

-1 






9 Z CJ° 


= -j 

ftt siz ~ 

ft 0^1 

- 


F/S 

=-* 

= 


= 

-i 




fto & , 


- <?%- 
P* c/y 

Rl^ 

r -J 

8h«i ~ 


=■ 

- PV" 
Px* 

««7 

=-i 

ffl* 7Z 9 


*** 

_ P 2 <^ 
Py * 

ftl t$8 

-3 

4,3 $ ~ 

ft&CI 


* Z 6T° 

■ 7T3^ 


= -l 

^ ~ 

RO^ 

*4 

— y .. 

-1 





~l 


■1 
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Cl 2 .si) 



1 2.4 


TRIANGULAR ELEMENT DO MATRIX 


[do 


E* 

i-V 2 - 

o 

o 

o 

o 

J - V 2 - 

fjt 

>- ys- 

o 

o 

o 

o 

\ 

o 

o 

aof-v) 

O 

o 

o 

o 

O 

O 



o 

/ic/-v*) 

n(/-v z ) 

V 


O 



EE Z 

O 

o 

o 

/ZC/-V 2 -) 

/Z.O-V 1 ) 


o 



O 


o 

o 

O 
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TRIANGULAR ELEMENT HMN EQUATIONS 


12.5 


Depending on the option exercised, the HMN equations have the form 



OH 0 

or 6> X <e 
or /8*l8 



OHO 
Or &X/£L 
or /QX./Z. 


OHO 

6x6 

Of / 8 x 6 


(1Z.-S3) 


where [Hi] is the product of a triply subscripted matrix and the column of total 
deformational rotations accumulated in the previous n-1 steps. 


Hi ~ Mlijk 0, 


4* Cn-0 




The [ H ] matrix must then be inverted, expanded and reordered .to form_the con 
straint equation. 


The C H ] [HQ] and Hi... 

ijk 

Terms that are not written 


matrices are presented on the following pages, 
explicitly are zero. 


REPRODUCIBILITY of the 

original page is poor 
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NUMBER 

JBUJFfA/iJ 

REV ITR 


AC„ 

U K- 

K‘) 


M,*' 


Z(V, 1 ? 

*r) 


IZ*>” 

eu. tn 


M m ) SA /1 


i(sf- 


-ey/'k *, 

. c*) . . 

5 i Sm o( l 


ZS, CoS A, 


~ 5 f L.ief, 


*■- 

V*) 


- tz*S 
>) ss? 

sLs.“>- 

3 .'*) 


5,‘0 95< a 

sa ( " 

3“') 


1 s/*) 


v ,j ) 


s/») sir 


fl* 4 « 

y<V a - 

m w ) 


-v") 


2 (V/ M - /Za 4 “- 

v/«) av ," 1 

rfa"’- 

^"V 


2 <y/'t 

W") 


111 t» oj 

izy, - -25 , j;, *, a 5 ,c«rf, 
85 ^' 


- 56 / 

s.‘”) 


< f»i , 


**Vf 


&(S, oi uad,- 


W/'Vnd.i 
- 5 i° Lu«t,* 

'W/'W, 


_ to 

■ 2J, J<‘n if, 

<« 1 . _ 

•a *n“i 


2 6 ( Vrt X,+ 
Af 4 '^ ot, 


2-Gi ewrf,- 
5 Z 


tzK, 1 "- ZiS^toX, 

«V/ n - ♦ /U/’W.) 

"Sj Co 1 »j ♦ “S?Sm<l s - 


Zis?S,\^ 

^ *1 (»,_j \ 


original page is poo. 




i c< 


NUMBER 
REV LTR 


4 <9w; 

A 


4 &&), 

4<&£ 

A0J-,, 

O.L, 


*r 

4r 

*c 

^ - 

ft, Fm 

ft, F-c in 

ft, ^*3/a 

ft ft it 

ft, FfyX 

ft, ft *iZ 







ft t ft at. 

k h*. 

fit F3ZZ 

ft Z. ft HZ Z 

Pt Fra 

ftzftii. 







ftj F <3* 

fts Fxsi 

kh,L 

R>F m 

rtfti 

ffjft^st. 






r 

fttFiz. 

» 

fit Fazr ! 

ft F jje. 

fl,r m 

ftFsn. 

ft, F i/z 

- 


'^Sii »t, 


"15, 5i'hS< ( 

3S%& ol, 

f 

1 







S,1;,«, 

”S, 

f) 

5, told, 

ftj,**. 

fltftvzz 

ft~i_F 3£i 


ftz.F: Sit 

ft-t- ’'iZl. 

~3S^*sihe! i ■ 




"5 £ 5 iVj o< z 

$^Zc6d- 


; 





*W £ 

~S*\oS ft 

• 


5 £ *kd t 

^AaSOtj. 

ft* *~i3Z. 

&F«t : 

^ jyi 

ftj F-iiz 

ftt ftrtc 

ft? F&j t 

.(3). 

”“j ?ih o(j 

Sj to5P<j 

"5 3^* 5 iiiodj 

3J^L>S0( 3 









-(?> 

S 3 j'y, 0(.j 

^Cs) 

Sj £osol s 

4 Singly 

.«) , 
"5j dc3«j 

‘ 




ft, ftz,F 

ftft 3 u* 

ft, ftltz* 

ft) 




r ft Fee 

+ft,ft?rc 

+ft, ftytx 


r ft,Fc-.t 



ftz ft IV- * 

ftzft i(F r 

ftzF^xT 

fti. ft J"tl + 

AftciS 


+ Rl ft IZT. 

^ ft l ftlx-Z 


+R t F<i zz 

* fti ftnz. 

^zft^z. 

MfV'ttti) 

-st^Cc-stt-h 



» « 


J __ 

, # 

*N z ^Su t d-z 

ft$ft UZ r 


/?3 ft?ZZ r 

ft's ft Hit 

ft's ft S'? l.^ 

ftfuS 

S } 'cesd-r - 

*\ ftit 



T ft 3 ^iSt 


*ft} 1 .•a.’ 

l lJ ’ * 

A'j ? n <A) 


<■ 1 ^ 

5 , • 


5 |jii 4 , ■* 


AV^d. ^“caa*, 

'■5/'ioi at, + -5 , W j,h od, - 
*A/| 6 ) SiVi /j/'isastrt, 


rtf) 

"S 3 < 2 c'J (* 3 * 

+>Vi m s;»u* 


•3£s*'s,nd t i- 

A/^WJ 

'ij l V(« jSj." 

A/j*’ V/1S 

„«} ✓ <•»)* , 

-S? SiViDljT J(5j 

t Mi dto o£j N/^sUii,) 

-n t/j ' "4 1-5 5<3 * 


■WjAfSld,) 


■* s(sS** r 

1 ‘ *<d*«0 

- •# “if s!» d. t ~ 

, ^‘dcaK, AJ/ J VcxS(rt, 

5 i'tc>^ L - S i tL fi n (r! i + 

*-ht t CcS*x. 

,l« , , - fc> , , 

"5 c Co6p( t 4 -5 4 

+ A4 ft Vh 3Ci A/l^'CC-S 9(1 


Of 4S02 81 


J/7f 


The HI matrix is written as a triply subscripted variable. 


ni.jk 1 

^*,2 4k ^ifcw Fjiz 

H 1 aj k. 

> I 

F P 

1 KIZ. 1 JIZ 

Hljjk 

ft *4 

F *22. h jn F^Z.\ FjZ2. 

RJhjK. r 

* 1 
f~ K 8.& H/ 2 £ 

hi I? jn. " 

< • * * 

F K3Z f" Jjfl ^ Fj 3& 

Hi <» j K 

* 1 

Hi 7 j fc 

F k>z Fj„ +Fk u F Jt z 


F f * 

‘K/2. 1 J |2 


Fja, f F*ai 

ffi'.jk = 

F KZ2- Fjzz 


F^sz Fj 3l + F h3 , Fj,z 

H2, I; k ~ 

F IC.3Z F j 32 

- 

Ftu . P j'// f F Wit fj/s f Fk., 2 p;„ + Fjo, Fj)t. 


K.Z F\^ hjhz 

H 1 /S' j k " 

i _ i — * “T * 

Fjzi f F kzt Fj zz 1 F,ai Fjzi i h*i i 


F K2t Fj*.* r j*-i 

j — . i 

httnjK 

F kJZ, Pj3| + 0** f F , <72 ^/3f + ^*3/ 

Hi * 

F *3£ Fj sz + F k3 z Fj 3 ?_ 


(,11 .vO 


The variables used to define these matrices are: 
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■RFPROflUCIBIUW °P ® E 

Sax, page is poor 


/? = ft tof- At''> ft (-wf) *M. t0 
ft = f>,(£ - s?) * hlrs*) ' ft 4’ 

ft - <**■”) 

ft --■• ft 5?' * ft (sf'-s , 1 * ) * ft 6 

ft - ft (wf ) * ft A//’ t ft ( A/, 01 ’ - ft 

ft .. R,c-sn*fks* 

See Table 2 for an illustration of the following area coordinate transformations. 





W? = ~n. 

Q (f) n 

Jg - ~ ^3 


nr 

~J 3 ~ 'dj 


wr = nr 

5T= -Jf 


0 

5*>r 

(12 30) 

< = "4' 

< 13) a fc) ff 

->j ~ •" -<4 , 


A/f'- -r?f* 

5, e>r 


.,(31 _ *) 

/ft - ft 

sf ■ -4^’ 


4 

■s~r 

If 

o 

if ~~ jf 
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• wr 

• 'i i - 9S ’m 

*m xk aevjfcte. 


\» fajtk 



- 

— T 5 t0 
e J 1 

^■331 

- 

— T 3 ^ 

£ 

‘SB/ 

= 

— J < tl> 
a J r -->2 

F nz. 

r 

. -L_t < cr; 

2 J l J ( 

F 332 

= 

1 -r- 

F's-zz 

- 

~tJs3z* 

^131 

- 

~T 3 ^ 

F sat 

— 


*V W 

- 

? < <l ^' 

^13Z 

- 

c ^ 

-y 3 s 3 - 

F SZ2. 

— 

~ p <r 

12 

— 


F 

1 izo 

- 

<<** 

"93 S 3 

P-3lO 

r- 

y,^ 93 

^jio 

- 

-rf 

£*/ 

- 

0 < <- 2 > 

^ 3)1 

- 

2 3 C, V- 

(~S-3I 

- 

(?) 

^ ' 


(iq.4®) 


ti i 


J r 5 


c/'y 


J—T < C1 ' ) 
a 'J I 0 3r 

I _ . a') 


lO 


a) 


O') 


«) _/_ 

/ + a 

- 0) o 

>, 

. CJ> 


* 

->a 

-03 , 


. 0 ) 


' iVL 

F~ 212 

^S3Z 
f Zll 

^j| 

fell 


a. Jj* J 3 


^4 


c/i ✓ CO o «r 0^ _J_ T c ^ 

2^3 


< u > O < 

■+ *^2 ^ A 


& 


zJx$3 


- -£-I*S 


(0 


0) 


= iVS 

r *T 5 t 


a) 


117 


^Z.1 1 

J ' T <r ^ 
i •Jz- ) i 

As- 

- 1 T 


1 _ „ (i) 

F &2Z 

“ z 5 Z 


. a) 

fit* = 



(Vi. ,4 c-) 

(cuHT.) 



ST 

"XA + ‘^Z2^Z 

F^jo 

= 

*3zS 3 ‘*-Xl.fF 

fzu 

- 

z* 3 S, ta * a *»“' 

II 

-- 

2 X ^3 5 3 + a J-v ^2 * 3 5, "2-Xz-)Sz 

F r 63t 

- 

( 3 ) L*} _, cs) , ) T < c3 ^ 

-i% 3t S 3 + 2 *^ 

f ZZt 

>-• 

c<A 1 - - . «■> . -C *- 1 

-Ex^S, ' 2 TzA *xA + #3 5* 


~ 

, fA ' -r x c b v <, ^ J, v cA 

~Z%^S 3 ' z J^I "'^ v 3^i 

f'&ZZ 

= 

< li’) / T <• < 3 > <- C ^-x < £3) 

2-x^5 z 2 , ^c. 5 3 

F Z3) 

— 



= 

-ax 3 ^ VT,5f 

Ku 

- 

_ ^ cb , 1 *r < tb 

2 - *32, ^ "* £ 

Fzzz 

- 


F ^iZZ. 

- 

2x 3 5- J ' i) -aJ^“'-X 3 5 ( <t> + X I 3 4 £l ' ' 

F <*iz 


- £ X 3? 5^ fi - a J t 5/° 4 Xsa 53^ ~ X* 5, 


(') 









The parameters F |ml and F.^ are defined to be the same as F. ml and F.^, 
except replace S.( m ) by N.( m )o 
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REPRODUCffiMTY OF '! 

,-vmritKrAT. PAGE IS ro> 


(V 


f * } 

/ &* 
i j 

: «• i 

f j 


- - 1 Q 

T * i i V 

'V 4 , f 

1 1 

i 9 
l vJ 




( 1 - 2 -4 5 *) 


where [T v 3 is the same as [T w ] except ^ v 66 ^ v 99 

These matrices [l u ] f [ T y 3 and [T w ] must be reordered to form 


instead of -1 . 


yl- yy 


IU-4G) 
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1 1.7 


QUADRILATERAL ELEMENT GO MATRIX 


G 0;,j = GO, J j 4 .|= , = GO;,j* 3 = G 

/=■*, 3,5,1,11: 

GO;, j =* 


= i h(i-T) 


J - )(. , 3C, 

n = V, 8 . 


;= 2,^6, /0,/2; 


GO;,] =£7nH„?)fM*27h?) f 


j= i,u ,zi,n 

n = 1,3, £■ ,7 


‘ 5 ?n o- r) 

= - 


f h 7 

l ^ - 2 ,6 ) 

C j- It >3t ? 
in* 4,6 / 


NOTE: j indicates column index in [GO..], 

n indicates value of coordinate at node number n. 


f = 1, 11, 21 .... 

n “ 1, 3 , 5 .... . 


etc. means use n — 1 if { 

" n = 3 if { 

" n = 5 if { 

etc. 


1 / 

= 11 


21 , 


(l2,4“f) 
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REPRODUCIBILITY OF THE 
ORIGINAL PAGE IS POOR 


[GO] Matrix - AZI Function Contributions, cont'd* 


GO; 0 = fiO.-. jtl 

i --7,8: 

-- i(<- f 

-- i (i-r)(uu?) 


J z /v.2% 3'i 
n -- t,3 4 sr, 7 w 


J : L 
n •= 2 


,2*0 
, <5 / 


r j- ^3<?? 

lr?* ^8 J 
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[GO] Matrix - HMN Function Contributions 

Let G = G(f,7) represent u 
G = G(7/^) represent v 


(l2.*Vu) 


G-I a n a„^,7) , S-- I j) « £ £„6„ Cf ,7) 

n*i n*i J rui 


^6 ^ ^6. 

M "ib, Qnl T * = 


£i„ 


<}£ 


r7 




P ' 




n G n 

i (f-f)(i~v) 

2 - ( f-f 3 )(/+y) 

3 

v (i-vxi-vn 

<r ZO-PW-'? 3 ) 

a i(P<){\-r)(T) 

7 i(h 0(i -?')(>?') 

a UyiKi-TK^) 

<* -4^wX/-?*yw 


G-n 

(l-T)OP) 

(7-? 3 ) (I*f) 
z(y-p)(!-f) 

(i-T)(i-p) 

z-O-TKi-y) 

i(?-0(i-P)(P 
i (?*/)//- r)(T) , 

nv-oo-riif) . 

u^oo-m) 


J? 
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[GO] Matrix 


HMN Function Contributions, cont'd. 


n } rn 

^Gr? 

hGfn 
& V 

1 

0-3f*)ti'7) 

(i-3y)(/-f) 

z 


(l-3? 2 )0+f) 

3 

Z( 1-30(1-7*) 

z(i-3T)o-r) 

V 

-zjd-T) 

-27(1 -fV 

S’ 

-*/ ft?- 

-Vi f - V) 

£> 

5 ^ 2 -!7 v J 

ur-y) 

7 

i (r-r; 

N 

1 

^ * 

8 

4 - 7 5 ) 

iif-n 

9 


ill- f) 


^ G r> 

<^G rn 

r\ } rn 


M 

1 

-if-V) 

-(7-7 3 ) 

z 

(f-D 

(v-r) 

3 

-v(j-y) 

-ifi?-r) 

V 

-z?(t ' t v 

-2)(J-V 7 ) 

S' 

2(l-37 z )(l-V) 

z(l-3p)(i-V) 


(7-zrxf-i) 

(0z0(7-i) 

7 

d-zrxf+i) 

(f-zyxv+o 

8 


j(i-30)(7-O 

c \ 

y fl-3 7*7(0! ) 

0 1-30(7+ )) 


Cl 2.51) 


(VZ.S2) 
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[GO] Matrix - HMN Function Contributions, cont'd. 
Columns 41 thru 8: 

j, m, n correspond according to following table. 


j - -y / 72 i3 77 75* 

n ~ V 5 3 9 

/77 s- " " " ^ 

j = 50 5/ 52 53 57 

n ~ Z 3 - - 7 


76 7 ? 75 77 

6 8 9 ~ 

55 56 57 55 

6 


tr? s 





GlO 




GOj'i - 




2 J ~ ^ 


^ ^ <^G(n 

G0 *j = rr 


GO 




m 


" a^ 


(±2. Si) 
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12.8 THE TRANSFORMATION J 
The transformation required is 


OF TH$ 




where the matrix [ J ] is known from Equation 


[J] = 


At 


Am 




At 


Az. 


a? 

c 

dx 



a) 


*5 


bu 




r 

Am 

As/ 

)W 

by 

a* 



*f 

Tf 


b* 



-l 

Au 

i/ 

^ (O’ 



4 

= w 

A'? 



[ A u. 


<J <U" 




Acs 

1 F2 

L_ 

^2 



AS 


J' 1 


Therefore evaluate [J] at each of the integration points, and invert 
to obtain the elements of [J ] at each of the integration points. 




Cvt.sO 


l 




[j]*' 


1 

■J.. 

» 

j IZ 

j>3 


"i v, x 7 " 


i 

Jzi 

Jzz 

Jm. 




> 

M 

4* 


X°-Jg* £ _ 


(l-l.S-c) 
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12.9 


QUADRILATERAL ELEMENT AO & A1 MATRICES 


J - ) V*. > j 4 Q > ^ ^<x, Qtj i j J S9 ^ 

\A£ \ - ~\ 4|a“'j 


AE; = ( RO ; j + FUjjk&t; )A0J 


#0 V = 

2 


ft If?? 

flOi? s 

dux 

ax 

0 


#0*;/ = 

2 



^2,6 = 

g>0> 

J s 

V 


#^2- ~ 

1 



#03,3 " 

1 


R1 J^flr 

RO 3<r z 

; +• 



& 

r' 

ii 

+ 

*or* 

Jx 


«v = 

■l 

/?4, = 

«kr° 
" 7T 


j 


dt*X° 

~ *V 

Ro^c. 1 

l 

» 

- 

~ ax 

O 

>j 

»t 

-i 

^r = 

-^r° 

a? 


(l/Z.fcl) 

2 

1 

('I'i.tl) 

1 

1 


pi w * ; 



Ri v r'^ 

fUsA*. ~1 

PUv" 1 


PI eft ~ 1 


:j 
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*5* 


R0 i.\ 

' ^x 8 - 


--i 

91 6, i, n 

ffo^ 

a^0J° 
J* c)y 

#!<.,% 3 

*-1.0 

91 (.M 

1 


fio b) ,r- i 


d z ur° 

HO 7fc ~ 

t *,H * '1 

9 1 7 p io ~ “ 1 

po^r-i 


fi 1 7 , C, '2 ~ 10 

Jlr° 

flOg ! C 

«)V° 

' “ a**- 

**w - 1 

91 s x tjiz ~ 1 

„ - $*•* 
90 83 “ ~ lcj x 

m tp? ~ ~i 

91t l \‘i = ’1 

- d'uT* 

R°t? " ~ t* i 9 
ffo V-'l 

~ 1 

91 g,i,n ' 1' 

fog^z-l 

91 g t , 0 / ' ’1 

9 1 $ ( i } to ~ ^ 


(12.. us) 


(17. . C ^ 
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1 Z „11 QUADRILATERAL ELEMENT HMN EQUATIONS 
HMN Equations LHS 



bj Q.-j Q ^ b ( b^. 

0 1 o 

1 i. j 

loo 

-i i -1 O 1 

-i i 1 -1 o - ■ 

10 -if ■ i - 1 

O ~1 ~"Zf 1 

0 1 o 

JL ± 1 

& 2- * 

1 0 o 




Then defining 

<fm7n j the values of ^ and ^ at node 


m 


X31 


% s ^ ^/,^7 

', 3*3 i 


2. 6 <SJ* 0* y«)()+ %) 12. U° £ (f*X )((* 7«) tf'2,6 

2^0 ****** ^ 




r? 


z±H<a:-z*or:ti 

w ■ 


t i cj° n - i rt u: t 

z,<. 4 V /7 


', 5 , 4 ? 7 


*S 3 ,s ;7 


? 2^ L L 0'7 n )(i'%)-t^^ f^O-XK/-^) w= i^j? 7 

'i 3 i^i' Z,(, ' 

f iu: (<-%)(!-?*) -i Ti*,: £o-v.)ib%) 

* i t ** t ' 


M=E ,6 


6 ^-h) 





tA-S-jQ 


n-->, w 


hA - 2 , &> 
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Ziz [i 1 1 + U 1 1(\> l 

*)]«•„ 

*- 1)0 

Z i[% (>*■&) *U (»tn)lvl 

i o *» *» 


TIIW0DT3CIBXUTY OF 'FEE 
ORIGIN AE PAGE IS POOR 






H Ja[jn 7n ?«U'U)* f* 0~ fn)] (*3°n 

I Z S' 1 

-t (<-%)* tX(i-U) 7^ 

7/§ 

H 8 

'I Tb[fn?n(l-fM) 4 ^7*(/-%)]v« 

/; J, 3^ 7 






21 m ^ 7*7«("%)6 f ti)-£&7M(i*k)( l *f'*)vZ * = '.W 

1, 3,^7 

E i (i*$J(/ f?„)cj° -i.~k7* {w$„)0<-jy)cj« m = %8 
w m 1 


Z zZ ?n ?» <UJ° ~ I 7t7;^ 

'/W 

5P _L 9 -LnZ c 

2 ~ & <^r„ - Z /2 7 ; &/„ 


M-'.W 


n~ v,B 


>>W 


Ol 27 ft 2 7 m (i-%)(>- - fjr iZo-fnjO'Dv: n-~ 1,3,57 

J 

/ I z(>-h)(i-UW„ -t7zt*(i-i)(i-U^ we 

v %s w> 
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HMN Equations RHS (HI) 


Note: H1,-H1 3 , Hl 4 -H1 5 , Hlg-HI,,, 

by changing sign of (1 tx.) or (1 tyi) 

i ' \ 


HI 




ZH fn frt ( l*Vn)( l*%) 

t(l* %)(>*%.) 
fk fn ( l r %)(l + '?n) 

4 


n,n * 1,3,5 7 

/7,M ' <2,6 

n- /j 3, S, 7 M= 5,0 

n - 2, 6 Af 1 /, 3, i~ 7 



n,rt= 

M r 2 ; £ 

n = / ; 3j5;7 r 2,6 

n 8 2, 6 M - l) 3 jSj 7 

M ~ I, 3,S~, 7 

/)j H = <2; ^ 

r? - Ij 3, 7 n - 1 ?■_, & 

r\-Z } (a r\ 7 1,3,$ J 

n ) M : 

- <2 j 6 

7 /n = 2,6 
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?« fa 0-%) +1 n \ ( l -%)] 

HI H^O-V* %(!-%)] 

\ rdv %( i -%)> fs?n ('-%)] 
fa (’- fa ) ■'$«'?« 0-7*)] 




A W*0*1n)(l+U) 
) &( i * fa )0* fa ) 

-fc ft ("■$»)(** fa ) 



- 




h,rt - 1,3,57 
n , m = H , 8 
n* 1,3,57 n*%8 

rr-%8 m-/ 4 3 y S' 7 


/ 7 jM= IjS^y 
Pjrt - H t 8 


I?* («W 


n -%8 ^ 1 , 5,5 7 

Ci'Z-^s) 


n , M r % 8 


t,<or<r.) 


rr-m<? A 1 =t 0 
/ 1=^3 


= /^7 
/ 7 jAl ' 9,3 

fA ~%& 

p*%s h ' l&fy 
n,M » /,J ,^7 

/ 4 ^,g 

n *%& M^ 3 ; 5; 7 
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REPRODUCIBILITY OF THU 
ORIGINAL PAGE IS 



* fjZjSJl 

Ciz.ca) 

/),M * *t,6 tco “ r -- ) 

? 7 * 
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